REPORT  NO.  NADC-82253-30 


S**8 

H  A  D  C 


IMPULSE  RESPONSE  OF  THE  OCEAN  FLOOR 
NONLINEAR  TECHNIQUES  FOR  MEASUREMENT  ENHANCEMENT 


Thomas  Gabriel  son 

Sensors  and  Avionics  Technology  Directorate 
NAVAL  AIR  DEVELOPMENT  CENTER 
Warminster,  Pennsylvania  18974 

April  1983 

Final  Report 
Task  Area  No. 

ZR  01108 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


DTIC  QUALITY  INSPECTED  3 

Prepared  by 

NAVAL  AIR  DEVELOPMENT  CENTER 
Department  of  the  Navy 
Warminster,  PA  18974 


NOTICES 


REPORT  NUMBERING  SYSTEM  —  The  numbering  of  technical  project  reports  issued  by  the  Naval 
Air  Development  Center  is  arranged  for  specific  identification  purposes.  Each  number  consists  of 
the  Center  acronym,  the  calendar  year  in  which  the  number  was  assigned,  the  sequence  number  of 
the  report  within  the  specific  calendar  year,  and  the  official  2-digit  correspondence  code  of  the 
Command  Office  or  the  Functional  Directorate  responsible  for  the  report.  For  example:  Report 
No.  N ADC-7801 5-20  indicates  the  fifteenth  Center  report  for  the  year  1978,  and  prepared  by  the 
Systems  Directorate.  The  numerical  codes  are  as  follows: 


CODE 


OFFICE  OR  DIRECTORATE 


00 

01 

02 

10 

20 

30 

40 

50 

60 

70 

80 


Commander,  Naval  Air  Development  Center 
Technical  Director,  Naval  Air  Development  Center 
Comptroller 

Directorate  Command  Projects 
Systems  Directorate 

Sensors  &  Avionics  Technology  Directorate 

Communication  &  Navigation  Technology  Directorate 

Software  Computer  Directorate 

Aircraft  &  Crew  Systems  Technology  Directorate 

Planning  Assessment  Resources 

Engineering  Support  Group 


PRODUCT  ENDORSEMENT  —  The  discussion  or  instructions  concerning  commercial  products 
herein  do  not  constitute  an  endorsement  by  the  Government  nor  do  they  convey  or  imply  the 
license  or  right  to  use  such  products. 


APPROVED  BY: 


DATE:  ^ 


/-( Pt'L. 


<r  ?. 


FREDERICK  D.  AMEEL 
DEPUTY  DIRECTOR,  SATD 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Data  Entered) 


REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

1.  REPORT  NUMBER  2.  GOVT  ACCESSION  NO. 

NADC-82253-30 

3.  RECIPIENT'S  CATALOG  NUMBER 

4.  TITLE  (and  Subtitle) 

Impulse  Response  of  the  Ocean 

Floor:  Nonlinear  Techniques  for  Measurement 
Enhancement 

5.  TYPE  OF  REPORT  &  PERIOD  COVERED 

Final 

6.  PERFORMING  ORG.  REPORT  NUMBER 

7.  AUTHORfsJ 

Thomas  Gabriel  son 

8.  CONTRACT  OR  GRANT  NUMBER(s) 

NAVAIR  DEV  CEN 

INDEPENDENT  RESEARCH  FUNDS 

9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Sensors  and  Avionics  Technology  Directorate 

Naval  Air  Development  Center 

Warminster,  PA  18974 

10.  PROGRAM  ELEMENT,  PROJECT,  TASK 

AREA  &  WORK  UNIT  NUMBERS 

Task  Area  No. 

ZR01108 

It.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

Naval  Air  Development  Center 

Warminster,  PA  18974 

12.  REPORT  DATE 

April  1983 

13.  NUM8ER  OF  PAGES 

55 

14.  MONITORING  AGENCY  NAME  ft  ADDRESS^//  different  from  Controlling  Office) 

15.  SECURITY  CLASS,  (of  this  report) 

Unclassified 

15a.  DECL  ASSI  F|  C  ATI  ON/ DOWN  GR  ADI  N  G 

SCHEDULE 

16.  DISTRIBUTION  STATEMENT  (of  this  Report) 

Approved  for  Public  Release;  Distribution  Unlimited 

17.  DISTRIBUTION  STATEMENT  (of  the  abstract  entered  in  Block  20,  if  different  from  Report ) 

18.  SUPPLEMENTARY  NOTES 

19.  KEY  WORDS  (Continue  on  reverse  side  if  necessary  and  identity  by  block  number) 

Impulse  Response 

Ocean  Floor  Measurements 

Nonlinear  Techniques 

20.  ABSTRACT  (Continue  on  reverse  side  if  necessary  and  identtfy  by  block  number) 

This  report  describes  a  procedure  for  enhancement  of  impulse  response 
measurements  of  the  ocean  floor.  When  these  measurements  are  made  with 
explosive  charges,  the  source  waveshape,  being  long  and  oscillatory,  obscures 
much  of  the  detail  in  the  received  signal.  The  effects  of  this  source  wave¬ 
form  can  be  eliminated  by  various  forms  of  deconvolution,  one  of  which  - 
homomorphic  deconvolution  -  is  very  flexible  in  its  application  and  does  not 

DD  ,JANM73  1  473  EDITION  OF  1  NOV  65  IS  OBSOLETE  UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Data  Entered) 


UNCLASSIFIED 

SECURITY  CLASSIFICATION  OF  THIS  PACE  (Whan  Data  Entarad) 


20.  Abstract  (Continued) 

require  foreknowledge  of  the  source  waveshape.  Homomorphic  deconvolution 
is  a  nonlinear  process  based  on  filtering  the  logarithmic  spectrum  of  a 
signal.  The  special  problems  such  as  phase  reconstruction  and  signal 
conditioning  are  discussed  in  the  context  of  ocean  floor  measurements 
with  explosives,  and  a  computer  code  is  presented  that  generates  and 
checks  the  log-spectrum. 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Data  Entered) 


REPORT  NO.  NADC-82253-30 


LIST  OF  FIGURES 

Figure  Page 

1.  Basic  task  of  nonlinear  signal  enhancement  .  4 

2.  Verification  of  deconvolution  process  . .  5 

3.  Complete  system  for  ocean  floor  impulse  response 

measurement  and  validation  .  6  ' 

4.  Fourier  transform  path  in  complex  z-plane  .  17 

5.  Construction  of  response  for  a  single  zero  close  to  the 

Fourier  transform  path  .  19 

6.  Magnitude  and  phase  changes  in  spectrum  for  transform  path 

that  passes  near  a  zero  in  the  z-plane .  20 

7.  Four  adjacent  spectrum  points  of  a  properly-sampled  real 

signal  plotted  in  the  complex  spectrum  plane  .  24 

8.  Actual  spectrum  progression  between  the  four  points 

in  figure  7 .  25 

9.  Sampling  of  phase  and  phase  derivative  for  a  typical 

real  signal  sampled  at  the  Nyquist  rate .  27 

10.  Various  curve-fits  of  spectrum  data:  (a)  straight  line, 

(b)  cubic  curve  fit  to  four  points,  (c)  cubic  spline  fit 
to  six  points,  (x)  cubic  spline  fit  to  entire  spectrum, 

(o)  actual  intermediate  spectral  points .  28 

11.  Example  of  homomorphic  deconvolution  of  a  received  signal 
(A*B)  made  up  of  a  source  excitation  (A)  and  a  multipath 

medium  response  (B) .  32 

12.  Reconstruction  of  explosive  source  waveform  by 

log-spectral  averaging  .  34 

A-l.  Configuration  of  Fourier  and  Laplace  inversion  paths 

in  complex  z-plane .  A-3 

A-2.  Configuration  of  Fourier  and  Laplace  inversion  paths 
in  complex  s-plane.  Shaded  area  and  numbered  points 

correspond  to  figure  A-l .  A-4 


JXTIC  QUALITY  INSPECTED  3 

i  i 


REPORT  NO.  N ADC-82253-30 


TABLE  OF  CONTENTS 

Title  Page 

List  of  Figures . ii 

Summary  .  1 

Introduction  .  2 

The  Problem .  7 

Deconvolution  Theory  . .  10 

Applied  Homomorphic  Deconvolution  .  22 

Signal  Conditioning  .  35 

Conclusions . 36 

References . 38 

Appendix  A  -  Transform  Relationships  .  A-l 

Appendix  B  -  Computer  Code . B-l 


iii 


REPORT  NO.  NADC-82253-30 
SUMMARY 


This  investigation  was  primarily  motivated  by  the  lack  of  good  measure¬ 
ments  for  acoustic  properties  of  the  ocean  floor  at  low  frequency  (below  500  Hz). 
It  was  also  an  opportunistic  study  in  that  a  large  quantity  of  raw  measurements 
of  ocean  floor  acoustic  reflectivity  have  already  been  made  by  NAVAIRDEVCEN 
in  most  of  the  major  ocean  regions  in  the  world.  While  these  measurements  are 
not  of  high  enough  quality  in  their  present  form  to  provide  any  more  than  a 
summary  measure  of  acoustic  reflectivity,  they  would  provide  more  detailed 
information  if  their  time  resolution  could  be  improved  significantly.  Since 
the  measurements  have  already  been  taken,  this  would  provide  a  very  economical 
means  of  acquiring  acoustic  property  estimates  for  various  ocean  floor  regions. 
The  major  objective  was,  then,  to  develop  a  technique  for  enhancement  of  ocean 
floor  reflectivity  measurements. 

Because  of  its  flexibility,  a  nonlinear  technique  known  as  hormomorphic 
deconvolution  was  selected  for  the  signal  enhancement  process.  The  measure-  ' 
ments  in  question  were  made  with  air-dropped  sonobuoys  and  small  explosive 
charges  and  the  explosion  waveform  has  a  prolonged,  oscillatory  tail  that  ob¬ 
scures  much  of  the  time-domain  information  in  the  received  signal.  If  it  was 
not  obscured  by  the  source  waveform,  this  received  signal  could  be  used  to 
isolate  the  individual  paths  of  propagation  that  are  probed  by  a  particular 
source- receiver  geometry.  By  means  of  a  nonlinear  transformation  basic  to 
the  enhancement  process,  the  convolution  of  the  source  waveform  and  the  ocean 
response  function  is  reduced  to  an  addition  and  then  the  source  component  is 
suppressed  by  linear  filtering. 

Two  techniques  for  enhancement  were  developed  during  the  course  of  this 
project.  The  first  -  a  complete  deconvolution  process  -  has  been  computer- 
coded  and  tested  on  real  signals  with  reasonable  success.  Although  noise  makes 
the  procedure  tend  toward  instability,  most  of  the  measurements  have  very  high 
signal-to-noise  ratios  (a  virtue  of  the  explosive  source)  and,  in  conjunction 
with  proper  preconditioning,  the  process  performs  satisfactorily.  The  second 
process,  based  on  averaging  spectra,  has  proven  successful  in  isolating  the 
source  waveshape.  This  technique  is  also  nonlinear  but  it  exploits  some  of  the 
properties  of  the  source  signal  itself. 

As  a  result  of  their  nonlinear  nature,  both  of  these  techniques  are  de¬ 
pendent  on  the  nature  of  the  signals  processed.  Although  the  concepts  developed 
are  generally  valid,  it  is  unlikely  that  these  processes  would  work  without 
modification  on  other  types  of  convolution  problems.  The  system  proposed  by 
this  investigation  will,  however,  be  thoroughly  tested  for  its  ability  to  en¬ 
hance  the  type  of  signals  available  in  the  NAVAIRDEVCEN  ocean  floor  reflectivity 
measurements.  This  "production"  work  will  be  done  on  NAVAIR  block  funding  since 
the  process  development  and  preliminary  testing,  funded  by  NAVAIR  DEV  CEN  Inde¬ 
pendent  Research,  has  been  completed. 
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An  ancillary  task  of  this  work  was  the  development  and  testing  of  a  math¬ 
ematical  model  and  accompanying  computer  code  for  the  acoustic  response  of  the 
ocean  floor  to  transient  signals.  This  task. is  fully  reported  in  the  second 
volume  of  this  report,  published  separately.  The  model  is  necessary  for 
interpretation  of  ocean  floor  reflection  signals  after  they  have  been  enhanced. 
Unless  the  individual  propagation  paths  can  be  isolated,  accurate  physical 
property  measurements  cannot  be  made.  Because  this  study  was  concerned  with 
low  frequency  transient  signals,  existing  acoustic  propagation  models  were 
generally  unsatisfactory. 


INTRODUCTION 


Since  World  War  II,  research  in  underwater  propagation  of  sound  has  pro¬ 
gressed  from  simple  empirical  studies  to  detailed  numerical  computations  of 
sound  fields  in  the  inhomogeneous  ocean.  While  this  research  and  the  attendant 
mathematical  models  have  spawned  sophisticated  means  for  predicting  acoustic 
propagation  through  ocean  water,  comparatively  little  progress  has  been  made  in 
modeling  propagation  through  the  ocean  floor.  In  light  of  the  importance  of 
both  advanced  sensor  design  and  fleet  ASW  operations  in  shallow  water,  at  low 
frequency  or  for  near-bottom  sensors,  the  role  of  the  ocean  floor  as  a  propa¬ 
gating  medium  (rather  than  as  a  simple  reflector)  must  be  understood. 

This  deficit  in  understanding  is  not,  in  general,  the  result  of  a  lack  of 
theory,  but  rather  the  consequence  of  a  paucity  of  reliable  measurements  of 
the  actual  acoustic  properties  of  the  ocean  floor.  Information  about  the 
acoustic  behavior  of  ocean  sediments  and  basement  rock  is  sparse  and  is  often 
based  on  laboratory  measurements  at  high  frequency  of  very  small,  disturbed 
samples  -  measurements  that  cannot  be  extrapolated  to  low  frequency  in-situ 
conditions.  Consequenty,  high  quality,  in-situ  measurements  of  ocean  floor 
acoustic  properties  are  needed  before  truly  representative  models  of  this 
boundary  can  be  constructed. 

Over  the  past  decade,  NAVAIRDEVCEN  has  made  measurements  of  sound  inter¬ 
action  with  the  ocean  floor  at  approximately  250  sites  around  the  world. 
Potentially,  these  data  can  be  used  to  answer  some  of  the  fundamental  questions 
about  the  influence  of  the  bottom  on  sound  transmission.  The  time  resolution 
of  these  measurements  is  not,  however,  sufficient  to  permit  direct  estimation 
of  in-bottom  properties.  This  resolution  is  greatly  reduced  by  the  oscillations 
of  the  gas  bubble  produced  by  the  acoustic  source  (an  explosion).  This  gas 
bubble  expands  and  contracts  rapidly  while  exchanging  energy  between  water-mass 
movement  and  gas  compression  until  an  equilibrium  is  reached.  Ideally,  a  very 
short  pulse  should  be  used  to  isolate  individual  energy  paths  through  the 
bottom  but  the  gas  bubble  oscillations  spread  the  explosion's  pulse  in  time, 
thereby  obscuring  the  individual  arrivals. 

Classically,  a  method  known  as  inverse  filtering  or  replica  deconvolution 
has  been  used  to  reduce  the  smearing  effects  of  the  bubble  oscillation.  In 
this  process,  the  exact  waveform  of  the  bubble  oscillation  must  be  known  so 
that  a  filter  can  be  constructed  that  reduces  this  waveform  to  a  narrow  spike. 
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In  this  investigation,  a  more  general  signal  enhancement  procedure  is 
developed  in  which  the  source  waveform  need  not  be  known.  This  method,  called 
homomorphic  deconvolution,  can  be  used  not  only  to  reduce  smearing  by  the 
source  waveform,  but,  in  another  mode  of  operation,  can  be  used  to  isolate  the 
source  waveform  itself.  Although  inverse  filtering  has  been  successfully 
applied  to  bubble  pulse  reduction, 2  the  homomorphic  procedure  is  inherently 
more  flexible. 

Essentially,  homomorphic  deconvolution  involves  transforming  the  time 
series  signal  into  a  domain  in  which  convolution  becomes  addition.  Once  in 
this  additive  domain,  the  bubble  pulse,  which  was  originally  convolved  with 
the  ocean  transfer  function,  can  be  suppressed  by  standard  linear  filtering 
procedures.  Transformation  back  into  the  time  domain  then  yields  the  desired 
ocean  response  as  if  the  source  has  been  a  single,  very  sharp  pulse. 

In  conjunction  with  this  signal  enhancement  process,  it  is  also  necessary 
to  develop  a  mathematical  model  of  the  interaction  of  transient  acoustic  signals 
with  the  ocean  floor.  The  model  must  be  sufficiently  general  to  describe  re¬ 
flection  and  refraction  within  layered,  inhomogeneous  media;  distortion  intro¬ 
duced  by  dispersion,  attenuation,  and  phase-shifting;  and  reflection  of  non- 
planar  wavefronts.  Without  such  a  model,  the  process  of  interpreting  the 
enhanced  signals  physically  would  be  haphazard.  Furthermore,  if  a  model  can 
be  developed  that  predicts  signal  propagation  similar  to  that  observed  in 
measurements,  these  measurements  can  be  associated  with  physical  phenomena  and 
quantitative  estimates  of  physical  properties  in  the  ocean  floor  can  be  made. 

The  interrelation  of  the  different  phases  of  this  project  are  shown  in 
figures  1  through  3.  The  fundamental  task,  that  of  deconvolving  a  measurement 
signal  to  generate  the  response  of  the  medium,  is  diagrammed  in  figure  1.  A 
by-product  of  the  homomorphic  technique  is  the  ability  to  isolate  the  source 
waveshape.  This  can  then  be  reconvolved  with  the  medium  response  as  shown  in 
figure  2  to  provide  a  measure  of  the  quality  of  the  deconvolution. 

Once  the  medium  response  has  been  measured  and  verified,  the  response  is 
translated  into  a  physical  model  of  the  medium.  This  is  the  desired  product 
of  this  research  but  it  too  must  be  verified.  For  this  reason,  the  transient 
response  model  was  developed.  The  model  of  the  medium  is  excited  by  the 
source  waveform  (by  computer  simulation)  and  the  resulting  theoretical  response 
can  be  compared  to  the  actual  received  waveform.  In  this  way  the  assumptions 
involved  in  building  the  physical  model  of  the  ocean  floor  can  be  evaluated. 
Thus,  the  procedure  as  a  whole  (figure  3)  has  several  feedback  paths  to  insure 
validity  of  the  results. 

This  report  summarizes  the  development  of  the  deconvolution  process  only 
since  that  was  the  most  time-consuming  task.  A  companion  volume*  describes 
the  design  and  operation  of  the  model  for  transient  response  of  the  ocean  floor 
and  also  some  of  the  results  obtained  by  this  phase  of  the  project. 
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Figure  1.  Basic  task  of  nonlinear  enhancement. 
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Figure  2.  Verification  of  deconvolution  process. 
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Figure  3.  Complete  system  for  ocean  floor  impulse  response  measurement  and  validation. 
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THE  PROBLEM 


To  better  understand  the  motivation  for  enhancement  of  seismic  measure¬ 
ments  by  explosives,  we  will  briefly  consider  the  process  of  measuring  relevant 
acoustic  parameters  of  the  ocean  floor.  The  most  important  of  these  quantities 
are  compressional  sound  speed  and  density  as  functions  of  depth  into  the  sedi¬ 
ment.  At  low  frequency  (less  than  100  Hz)  or  in  shallow  water  (a  few  wavelengths 
or  less),  the  shear  sound  speed  and  attenuation  may  also  become  important. 

The  depth  to  which  these  quantities  must  be  known  depends  on  the  geometry 
of  the  process  to  be  modeled,  the  rate  at  which  compressional  or  shear  energy 
is  dissipated  in  the  sediment  (the  two  attenuation  parameters),  and  the  acoustic 
frequency.  Since  viscous  attenuation  (constant  loss  per  cycle)  is  generally 
assumed,  the  range  of  frequency  to  be  covered  by  the  model  will  set  some  limit 
on  the  depth  to  which  we  will  need  to  know  the  properties. 

Any  experiment  designed  to  measure  these  parameters  must  take  into  account 
the  actual  process  to  be  modeled.  If,  for  example,  we  are  interested  in  low 
frequency  detection  of  a  single  frequency  source  in  deep  water,  we  should  try 
to  make  the  measurement  at  a  similar  frequency  and  in  deep  water.  In  particular, 
the  scale  (in  this  case,  the  dimension-to-wavelength  ratio)  of  the  experiment 
should  be  close  to  the  scale  of  the  application.  Attenuation  in  sediments  is 
not  linearly  related  to  frequency  over  wide  frequency  ranges  (as  viscous  attenua¬ 
tion  theory  applies)  so  the  measurement  frequency  should  be  as  close  as  possible 
to  the  application  frequency  for  attenuation  measurements.  The  experiment 
itself  should  not  disturb  the  medium  significantly  and  it  should  provide  suf¬ 
ficient  resolution  in  space  so  that  the  parameters  can  be  determined  as  functions 
of  depth  and  range. 

In  short,  the  measurement  should  be  similar  in,  at  least,  scale  and,  perhaps, 
fequency  to  the  system  to  be  modeled  and  should  probe  the  bottom  sufficiently 
in  depth  and  range  without  disturbance.  Under  these  conditions,  we  must  select 
a  measurement  method  that  will  provide  realistic  estimates  of  the  properties 
important  for  modeling  underwater  acoustic  sensor  performance. 

Historically,  sound  speed  and  attenuation  measurements  have  been  made  on 
thousands  of  samples  of  ocean  sediment  obtained  by  coring.  For  our  purposes, 
these  measurements  are  seriously  deficient  on  at  least  two  counts:  first,  the 
coring  process  produces  substantial  deformation  of  the  sample  thereby  changing 
the  physical  properties;  second,  the  measurements  are  generally  made  at  several 
kilohertz  so  that  the  attenuation  values  are  not  applicable  to  low  frequency 
investigations.  Relatively  little  work  has  been  done  on  scaled  experiments 
and  it  would  be  difficult  to  construct  a  scale  sedimentary  column.  In  any  case, 
attenuation  values  could  not  be  measured  reliably  since  scaling  length  downward 
means  frequency  must  be  increased  (in  order  to  preserve  the  length-to-wavelength 
ratio) . 
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There  are  two  major  types  of  measurements  made  in  the  same  environment 
as  the  operating  system  will  be  a  part  of.  These  are  continuous  wave  (CW) 
and  wideband,  short-duration  "shot"  measurements.  Although  there  is  current 
interest  in  detection  of  broadband  energy  in  the  ocean,  most  ASW  systems  are 
designed  for  discrete  frequency  emissions  from  the  target.  In  this  sense, 
the  CW  measurements  ought  to  be  more  realistic;  however,  the  CW  technique  does 
have  some  disadvantages.  As  discussed  below,  the  shot  method  can  be  consider¬ 
ably  cheaper  and  faster  than  CW  and,  furthermore,  the  CW  signal  cannot  be  used 
to  isolate  different  propagation  paths  (without  beamforming).  Therefore,  the 
received  CW  signal  is  a  single  sine  wave  completely  characterized  by  an  ampli¬ 
tude  and  a  phase  and  can  only  be  used  to  measure  the  total  propagation  loss 
between  a  source  and  a  receiver.  If  the  water-borne  path  dominates,  then  very 
little  information  can  be  deduced  about  bottom  properties. 

One  of  the  primary  reasons  for  adopting  the  shot  method  for  these  ocean 
floor  studies  was  that  the  entire  experiment  can  be  done  from  an  airplane 
using  sonobuoys  and  small  explosive  charges.  In  comparison  with  ship-based 
experiments,  this  technique  is  cheap  and  capable  of  surveying  large  areas 
quickly.  The  explosives  -  standard  MK81  SUS  -  are  reliable,  repeatable,  and 
result  in  a  high  signal-to-noise  ratio  at  the  receiver.  Also,  the  bandwidth 
is  high  (roughly  50-10,000  Hz  for  a  2  lb  charge  at  800  ft)  so  that  the  potential 
time  resolution,  after  enhancement,  is  good. 

On  the  other  hand,  precise  placement  of  the  source  and  receiver  is  diffi¬ 
cult,  the  source  level  at  very  low  frequency  (below  50  Hz)  has  not  been 
accurately  determined,  the  oscillating  gas  bubble  makes  enhancement  difficult 
and,  because  of  the  large  amount  of  energy  released,  the  initial  propagation 
of  the  pressure  pulse  is  definitely  nonlinear.  This  last  problem  can,  in  some 
instances,  make  application  of  the  results  to  linear,  CW  propagation  question¬ 
able. 


In  this  report,  we  will  consider  only  the  enhancement  of  the  received 
signal  by  removing  the  smearing  effects  of  the  explosion's  gas  bubble  oscilla¬ 
tion.  This  bubble  is  formed  under  very  high  pressure  by  the  explosive's 
combustion  products  and  it  expands  until  the  gas/water  momentum  is  checked  by 
the  surrounding  fluid  pressure.  The  bubble  has,  at  this  point,  over-expanded, 
so  it  is  forced  back  in  towards  its  center.  Since  the  gases  are  quite  elastic 
and  the  flow  of  the  surrounding  water  is  spherically  symmetric,  this  oscillation 
continues  for  a  number  of  periods  until  the  energy  is  essentially  dissipated 
and  an  equilibrium  is  reached  between  the  gas  pressure  and  the  static  fluid 
pressure.  Typically,  four  to  six  cycles  are  visible  on  the  received  waveform 
of  such  an  explosion.  Without  eliminating  these  oscillations,  only  gross 
properties  of  the  bottom  can  be  studied.  If  the  oscillations  can  be  removed, 
the  time  resolution  of  the  received  signal  can  be  greatly  improved.  Consequently, 
property  estimates  along  individual  paths  through  the  ocean  floor  can  be  made. 

Before  we  discuss  the  process  by  which  these  received  signals  are  enhanced, 
let  us  examine  the  formation  of  the  received  signal.  While  several  points  of 
view  are  possible,  we  will  treat  the  problem  as  a  system  problem  with  some 
excitation  (the  explosion  and  bubble  oscillation)  and  some  response  (the 
received  waveform).  The  system  is  the  ocean  and  its  boundaries.  One  form  of 
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excitation  which  we  will  not  consider  is  a  single-frequency  sine  wave.  If  we 
were  to  measure  the  output  of  the  system  for  enough  of  these  sine  waves,  each 
at  a  different  frequency,  we  could  completely  determine  the  system  response. 

This  would  be  equivalent  to  measuring  the  received  level  from  a  sonobuoy  for 
each  setting  of  a  variable-frequency  CW  source.  In  this  manner,  we  could 
measure  the  frequency  response  of  the  system  and  then,  given  the  spectrum  of 
any  arbitrary  excitation,  compute  the  response. 

While  the  preceding  method  seems  reasonable,  it  would  in  fact  be  difficult 
to  do  because  of  the  need  to  sweep  the  source  through  a  broad  range  of  frequencies 
pausing  long  enough  at  each  to  achieve  steady  state,  all  the  while  maintaining 
the  same  experimental  geometry.  We  can,  on  the  other  hand,  probe  the  system  at 
all  frequencies  in  the  band  simultaneously  if  we  choose  an  excitation  waveform 
whose  spectrum  is  uniform  (that  is,  flat)  across  the  band.  The  impulse  is 
just  such  a  waveform.  Theoretically,  the  impulse  has  zero  width  and  infinite 
amplitude  so  that  the  area  enclosed  is  finite;  however,  for  practical  pur¬ 
poses,  we  will  consider  an  approximate  impulse  that  is  very  short  duration  and  . 
very  high  amplitude.  With  this  impulse  excitation,  the  received  signal  will 
have  some  characteristic  shape  called  the  impulse  response:  a  fundamental  system 
property. 

To  use  this  impulse  response,  we  sample  the  excitation  and  then  replace 
the  sample  points  by  weighted  impulses  so  that  the  net  area  under  the  curve 
is  conserved.  As  long  as  the  system  is  linear,  the  response  is  just  the  sum 
of  the  impulse  responses  weighted  and  shifted  in  accordance  with  the  original 
impulse  sum.  This  weighting  and  shifting  operation  is  called  convolution 
and,  in  this  instance,  the  time-domain  waveform  at  the  receiver  is  the  convolu¬ 
tion  of  the  system's  impulse  response  with  the  arbitrary  excitation  waveform. 

This  convolution  process  occurs  any  time  some  physical  system  is  excited. 

Each  little  "bit"  (interval  of  time)  of  the  excitation  function  acts  as  an 
impulse  and  the  system  responds  accordingly.  The  output  is  the  resultant  of 
all  of  the  impulse  responses  delayed  in  time  to  correspond  with  their  respective 
bits  of  excitation. 

This  superposition  of  responses  depends  on  the  linearity  of  the  system. 

In  most  cases,  the  systems  with  which  we  will  work  are  linear,  but  remember 
that  the  explosion  generates  a  large  amplitude  disturbance  and  until  the 
amplitude  decays  sufficiently  with  range  the  propagation  is  nonlinear.  If, 
for  example,  the  explosion  wavefront  reflects  from  the  ocean  surface  well 
before  it  becomes  linear,  the  results  may  be  quite  different  from  those  observed 
in  a  linear  reflection.  At  this  stage,  we  will  have  to  be  content  with  the 
linear  theory  but  aware  that  the  explosion  may  introduce  errors  into  this 
assumption.  Also,  we  should  note  that,  while  the  signal  enhancement  process 
described  herein  is  a  nonlinear  process,  it  does  assume  that  the  system 
on  which  it  operates  is  linear. 

In  the  next  section,  we  will  lay  the  theoretical  foundation  for  a  decon¬ 
volution  process  whereby  the  impulse  response  of  the  ocean  system  is  extracted 
from  the  received  signal.  As  we  have  seen,  the  received  signal  is  formed  by 
convolving  the  impulse  response  with  the  excitation  waveform  which,  in  our 
experiments,  is  an  explosion's  pressure  pulse  and  bubble  oscillation.  If  this 
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convolution  can  be  undone,  the  impulse  response  and  the  source  waveshape  may 
be  separated.  Since  the  impulse  response  is  the  response  of  the  ocean  and 
ocean  floor  to  an  infinitely  short  pulse,  we  then  would  have  a  description 
of  the  propagation  with  potential  for  resolving  individual  propagation  paths 
without  the  gross  smearing  caused  by  the  gas  bubble  oscillations. 


DECONVOLUTION  THEORY 


Convolution  is  a  common  phenomenon  in  physical  problems  in  that  a  convolu¬ 
tion  occurs  any  time  an  external  force  is  applied  to  a  system.  Real  systems 
generally  have  some  lag  or  delay  in  responding  to  an  excitation  and,  perhaps, 
a  resonant  "ringing"  or  a  decay  that  prolongs  the  response  after  removal  of 
the  excitation.  Since  this  behavior  is  common,  we  will  review  some  of  the 
mathematical  tools  used  to  describe  convolution  before  considering  the  methods 
for  deconvolution. 

As  described  in  the  previous  section,  the  fundamental  excitation  is  the 
impulse  function  and  the  fundametal  system  response  is  the  impulse  response. 

We  will  use  the  following  notation  to  describe  the  excitation  and  response, 

6(t)  +  h(t). 

The  left-hand  side  is  the  excitation  which,  in  this  case,  is  the  impulse  (or 
delta)  function.  The  right-hand  side  is  the  response  (the  impulse  response 
in  this  example).  The  impulse  function  is  defined  as, 

8 (t)  =  1/d t  (dt  o)  if  t  =  o 

=o  if  t  f  o 

so  that  the  following  is  true, 

r~ 

J  5(t-tQ)  f  (t)  dt  =  f  (tQ)  (t1  <  tQ  <  t2)  .  (1) 

Thus,  the  impulse  samples  a  function  under  integration. 

A  shift  in  time  does  not  affect  the  system  response,  therefore, 

6(t-x)  h(t-x) 

and,  since  the  system  is  linear,  we  can  form  a  superposition  of  these  impulses 
as  follows, 


03  CO 

^2  fn  S(t-nAt)  -»•  J2  fn  h(t-nAt). 
n=o  n=o 
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For  a  continuous  weighting  function,  the  summation  becomes  an  integration. 


f(x)  6(t-x)dx 


cx 

f 


f (x)h(t-x)di 


but,  because  of  the  sampling  property  of  the  impulse,  the  left-hand  side 
reduces  so  that 

co 

f(t)  ->  J  f  (x)h(t-x)dx.  (2) 

o 

This  equation  expresses  the  system  response  mathematically.  For  an  excitation 
f(t),  the  system  response  is  the  convolution  of  that  excitation  with  the 
system  impulse  response.  If  we  know  the  impulse  response,  we  can  then  calculate 
the  system's  response  to  any  arbitrary  excitation.  (We  have  assumed  above 
that  the  excitation  starts  no  earlier  than  t=0.  This  is  not  necessary  -  just 
convenient. ) 


From  this  point,  we  will  leave  the  time  domain  and  concentrate  on  the 
frequency  domain.  Most  of  the  literature  that  describes  homomorphic  decon¬ 
volution,  the  particular  variety  of  deconvoluation  we  are  considering,  uses 
the  z-transform  to  establish  the  method's  validity. ^ > 4 , 5  s-jnce  we  are 
ultimately  going  to  use  the  fast  Fourier  transform  (FFT)  to  convert  from  one 
domain  to  another,  the  discussion  to  follow  will  not  explicitly  mention  the 
z-transform.  Instead,  the  Fourier  transform  will  be  used  and  connection  between 
the  Fourier  transform  and  the  z-transform  has  been  relegated  to  appendix  A. 
Hopefully,  this  will  make  the  process  description  more  palatable  to  the  non¬ 
specialist. 

Take  the  Fourier  transform  of  the  convolution  integral,  equation  (2), 

oo 

G(w)  =  jj'  f(x)h(t-x)e  ^dtdx 

’o 

Next,  expand  the  exponential  factor, 

oo 

G(w)  =  JJ'  f(x)e  lw^h(t— t)e  J  ^  T^dtdx 
o 

and  make  change  of  variables  as  follows, 

t'  =  t  -  x 
dt '  =  dt 

so  that, 

oo 

G(w)  =  f  f(x)e  lwtdxQ 
o 

=  F(w)  H(u) 

where  F(a))  =  Fourier  transform  of  f(t) 

H(w)  =  Fourier  transform  of  h(t). 


oo 

f  h(t ' )e"1ut'dt 
o 
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Hence,  a  convolution  of  two  time  functions  is  equivalent  to  the  product  of 
their  spectra.  In  symbolic  notation,  this  property  can  be  stated  as  follows, 

F[f*h]  =  F[f ( t ) ]  F[h ( t ) ] 

=  F((o)  H(tu)  (3) 

where  F  [  ]  =  Fourier  transform 

*  =  convolution 

As  a  consequence  of  equation  (3),  we  can  now  compute  the  system  output 
(in  the  frequency  domain)  by  multiplying  the  spectrum  of  the  excitation  by  the 
spectrum  of  the  impulse  response.  In  addition,  we  could  deconvolve  by  division: 
if  we  want  to  extract  the  impulse  response  of  the  system,  we  can  divide  the 
spectrum  of  the  system  output  by  the  spectrum  of  the  excitation  and  the  result 
will  be  the  spectrum  of  the  impulse  response. 

The  success  of  this  type  of  deconvolution  depends,  however,  on  two  factors. 
First,  the  excitation  spectrum  must  be  accurately  known,  and  in  these  airborne 
seismic  measurements  the  excitation  is  rarely  directly  measurable.  Second,  the 
excitation  spectrum  can  never  be  zero  (or,  for  that  matter,  small)  at  any 
frequency  in  the  band  of  interest  since  the  resultant  response  would  be  un¬ 
defined  (or  extremely  large)  at  these  frequencies.  This  problem  is  usually 
circumvented  by  adding  some  noise  to  the  excitation  to  fill  in  spectral  nulls. 

This  deconvolution  by  division  is  known  as  inverse  filtering  and  much  work 
has  been  done  in  this  area.-  If  the  source  waveshape  is  well  known,  this  is 
probably  one  of  the  best  methods  for  computing  impulse  response;  however,  we 
will  not  say  much  more  about  inverse  filtering. 

Homomorphic  deconvolution  promises  to  be  considerably  more  flexible, 
particularly  in  the  case  of  a  very  complicated  system  such  as' the  ocean  and 
ocean  floor.  As  we  will  see,  it  should  be  possible  to  isolate  not  only  the 
impulse  response  but  also  certain  parts  of  the  impulse  response  or  even  the 
excitation  itself.  It  is  this  processing  flexibility  that  makes  the  homomorphic 
technique  attractive. 

Linear  filtering  of  signals  is  a  very  useful  and  flexible  procedure  for 
noise  removal  or  signal  shaping  because  the  filtering  can  be  done  in  the 
frequency  domain.  If  we  consider  a  signal  to  be  a  superposition  (that  is,  a 
summation)  of  sine  waves  of  different  frequencies  and  amplitudes,  these  com¬ 
ponent  waves  can  be  selectively  removed  or  amplified  in  the  frequency  domain. 

The  basis  for  homomorphic  deconvolution  is  the  location  of  a  domain  in  which 
convolution  becomes  addition.  Once  this  domain  is  found,  the  highly-developed 
techniques  of  linear  filtering  can  be  applied  to  separate  or  modify  the  various 
components.  Thus,  homomorphic  deconvolution  involves  a  nonlinear  transformation 
of  the  convolved  signal  to  reduce  the  convolutions  to  summations  and  subsequent 
linear  filtering  to  extract  the  desired  information. 

We  have  already  reduced  convolution  to  multiplication  by  means  of  the 
Fourier  transform.  If  we  can  reduce  the  multiplication  to  addition,  we  will 
in  principle,  have  developed  the  required  process.  The  logarithm  function 


12 


REPORT  NO.  NADC-82253-30 


reduces  multiplication  to  addition  but  the  numbers  on  which  we  must  operate 
(the  spectrum  values)  are  complex.  Hence,  we  must  define  a  complex  logarithm 
such  that, 

Ln  (XY)  =  Log(X)  +  Log(Y) 

where 

X,Y  =  complex  numbers. 

This  condition  will  be  satisfied  if, 

Ln  (X)  =  ln/X/  +  i  arg(X)  (4) 


or,  if 

X  =  xe1'* 

then, 

Ln  (X)  =  ln(x)  +  i<[> 

where 

In  =  natural  logarithm 

arg  =  argument  function 
//  =  absolute  value. 

Another  equivalent  definition  is  possible  by  generalizing  the  definition 
of  the  natural  logarithm. 

In  (x)  =  /  <£ 

where  x  is  real.  If  we  substitute  a  complex  z  for  x  and  change  the  integration 
to  a  contour  integration,  we  have, 

Ln  (Z)  =  xfz^  (6) 

The  integral  form  of  equation  (6)  is  not  a  very  practical  definition  but 
it  serves  to  illustrate  a  problem  in  determining  the  imaginary  part  of  the 
complex  logarithm.  According  to  equation  (4),  the  imaginary  part  is  the  phase 
of  the  complex  number;  however,  in  the  exponential  notation  of  equation  (5) 
this  phase  is  ambiguous.  Any  integer  multiple  of  2v  can  be  added  to  the  phase 
factor  <t  in  equation  (5)  without  changing  the  value  of  the  complex  number. 

This  addition  will  change  the  logarithm  though.  Equation  (6)  is  a  contour 
integral  and  the  integral  has  only  one  pole  at  z=0.  Therefore,  any  contour 
will  produce  the  same  answer  as  long  as  the  contours  all  encircle  the  origin 

in  the  same  way.  The  addition  of  multiples  of  2tt  to  the  phase  of  the  complex 

number  is  equivalent  to  circling  the  origin  that  same  number  of  times  and  the 

resultant  value  is  a  function  of  the  number  of  origin  "orbits". 
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While  we  are  generally  accustomed  to  reducing  the  phase  of  a  complex 
number  to  its  principal  value,  we  must  not  do  so  if  we  intend  to  use  the 
complex  logarithm.  If  the  generating  function  for  z  is  a  continuous  function, 
the  phase  relationship  of  nearby  points  can  be  determined  by  examining  many 
closely-spaced  points  in  between  and,  in  this  manner,  the  absolute  (i.e.,  non¬ 
principal  value)  phase  can  be  established. 

Our  problem  is  complicated  by  the  fact  that  we  operate  on  sampled  signals 
and  sampled  output  from  FFT  routines.  In  this  case,  the  phase  can  change  so. 
much  from  one  point  to  the  next  that  the  absolute  phase  is  lost.  One  might 
naively  believe  that  this  could  not  happen  since  the  signals  are  properly 
sampled  to  prevent  aliasing  but  this  sampling  criterion  (the  Nyquist  rate)  is 
based  on  linear  theory  and  the  phase  is  a  nonlinear  function  of  the  spectrum 
value.  (The  phase  is  the  arc-tangent  of  the  ratio  of  the  imaginery  part  to 
the  real  part.)  While  the  sampling  may  be  adequate  to  define  linear  functions 
of  the  signal,  there  is  no  reason  to  expect  that  nonlinear  functions  are 
adequately  sampled  also.  In  fact,  the  phase  is  frequently  severely  under¬ 
sampled  in  the  signals  considered  by  this  study.  The  aliasing  that  results* 
from  this  under-sampling  is  discussed  in  terms  of  the  z-transform  by  Ulryclv 
and  Stoffa.5 

1 

This  then  is  the  principal  problem  with  homomorphic  deconvolution.  If 
the  absolute  phase  can  be  correctly  established  at  every  one  of  the  sample 
points  of  the  spectrum,  the  deconvolution  process  can  be  continued.  We  will 
defer  the  discussion  of  how  this  can  be  done  until  the  next  section  and 
continue  with  the  theoretical  development  as  if  the  absolute  phase  problem 
has  already  been  solved. 

Having  successfully  calculated  the  complex  logarithm  of  each  point  in 
the  spectrum,  we  have  effected  a  transformation  into  the  log-spectral  domain. 

In  this  domain,  the  respective  components  of  signals  that  were  convolved  in 
the  time  domain  are  summed.  Since  the  combination  of  components  is  now  linear, 
we  can  draw  upon  an  extensive  set  of  linear  filtering  techniques.  Usually 
these  techniques  are  applied  to  time  signals  by  first  transforming  the  signals 
into  the  frequency  domain  where  the  additive  components  (the  constituent  sine 
waves)  become  separated  along  the  frequency  scale.  In  this  same  way,  we  will 
transform  the  log-spectrum  into  another  domain  in  which  the  additive  components 
(which  were  convolutional  components  in  the  time  domain)  become  separated. 

Application  of  the  Fourier  transform  to  the  log-spectrum  generates  a 
function  called  the  complex  cepstrum.  Originally,  the  real  logarithm  of  the 
magnitude  spectrum  was  used  to  generate  the  cepstrum  thereby  losing  the  phase 
information.  This  was  still  a  useful  device  for  some  applications  but,  without 
the  phase  information,  deconvolution  can  only  be  done  in^some  very  special 
cases  (which  are  discussed  below).  With  this  early  work6  came  a  lexicon  of 
new  terms  describing  the  various  attributes  of  the  cepstrum;  however,  most 
of  these  words  serve  only  to  confuse  the  real  workings  of  the  process  so  we 
will  retain  only  the  word  cepstrum. 

Actually,  the  important  domain  is  the  log-spectral  domain  because  it  is 
here  that  the  convoluted  time  functions  first  become  linear.  The  cepstrum  is 
valuable  primarily  as  an  aid  to  filtering;  in  fact,  the  filtering,  once  developed, 
can  be  applied  directly  to  the  log-spectrum. 
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At  this  point,  let  us  consider  an  example  in  order  to  clarify  the  rela¬ 
tionship  between  the  log-spectrum,  the  cepstrum,  and  the  time  function.  The 
bubble  oscillations  of  an  underwater  explosion  affect  the  spectrum  by  intro¬ 
ducing  a  definite  oscillation  (periodic  along  the  frequency  scale)  of  the 
spectrum  magnitude.  A  Fourier  transform  of  this  magnitude  spectrum  (or  the 
log-spectrum)  should  show  a  discrete  peak  at  a  time  corresponding  to  the 
rapidity  of  the  oscillation  in  the  frequency  domain.  This  peak  does,  in  fact, 
appear  in  the  cepstrum  and  the  time  (the  scale  of  the  cepstrum  is  time)  cor¬ 
responds  to  the  period  of  the  bubble  oscillation.  In  general,  the  cepstrum 
scale  is  related  to  the  time  between  repetitions  of  an  event  in  the  time  series 
(such  as  echoes  or  multiple  paths  of  propagation). 

Another  more  specific  example  is  helpful  at  this  state.  Consider  a  system 
whose  response  is  a  series  of  echoes;  the  impulse  response  for  such  a  system 
is  a  sequence  of  delta  functions, 

m 

h(t)  =  X  a.fi(t-tj ) .  (7) 

i=o 

The  response  for  an  arbitrary  excitation  f(t)  is  then  a  sequence  of  weighted 
and  delayed  replicas  of  f(t), 

g(t)  =  f*h 
m 

=  X  a./f (t-r)6(r-ti )dr 
i=o  1 


m 

=  La.f(t-t.) 
i  =o 

where  g(t)  =  system  output. 

Furthermore,  the  frequency  response,  equation  (3),  is, 


/<»  ill 

X  6(t-t - )e-1 


or 


-Co  m 

i 

o  i=o 
m 

x 

i=o 

and  the  log-spectrum  is, 


-  i  cot 


dt 


G(u>)  =  F(u)  X  a.e' 


1  0)ti 


Log  [ G ( oo )  ]  —  Log  [ F ( oo )  ]  +  Log 


r  m 


X  aie-iwti 


1  =o 


More  specifically,  consider  one  arrival  and  a  phase-inverted  echo  delayed 
by  2t  seconds  so  that. 


m 

X  a-je-^ti  =  ]_  _ 


-i  o)2t 


i  =o 


=  2e"i(“T-^  sin  («t) 
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and , 


Log  [G ( a) ) ]  -  Log  [F(co)]  +  In  2 

-i  (cox-^-)  +  In  [sin  (cot)]. 


(8) 


The  last  two  terms  on  the  right-hand  side  deserve  some  comment.  The  last  term 
is  periodic  in  u  and  will,  therefore,  lead  to  a  peak  in  the  cepstrum  at  t 
seconds  which  is  proportional  to  the  echo  time  or,  in  other  words,  the  repeti¬ 
tion  rate  of  the  basic  excitation  waveform.  The  next  to  last  term  in  equation 
(8)  is  linear  in  w  and  can,  therefore,  lead  to  very  large  values  in  comparison 
with  the  last  term.  This  could  easily  swamp  the  more  important  final  terms. 
Fortunately,  this  term  can  be  eliminated  merely  by  shifting  the  waveform  g(t) 
in  time  (by  x  seconds  in  this  case)  so  that  the  two  impulses  of  the  h(t)  func¬ 
tion  are  centered  around  t=0.  This  term,  the  linear  phase  ramp  term,  is 
generally  present  in  the  log-spectrum  but  it  can  always  be  eliminated  by 
shifting  the  time  series. 

There  is  another  aspect  of  the  complex  logarithm  that  causes  a  consider¬ 
able  amount  of  trouble  and  that  is  its  value  for  very  small  arguments.  In 
real  signals,  we  are  normally  not  troubled  by  poles  in  the  spectrum  but  nulls 
or  regions  of  low  level  are  common.  Through  the  logarithm  operation,  zero 
points  in  the  spectrum  become  poles.  Thus,  any  regions  of  very  low  level  in 
the  spectrum  will  be  troublesome  in  the  course  of  tranformation  to  the  log- 
spectrum. 

To  visualize  the  relationship  of  zeros  in  the  spectrum  to  the  transform 
process,  we  will  generalize  the  Fourier  transform.  (Appendix  A  describes  some 
of  the  relationships  between  the  Fourier  transform  and  several  other  transforms.) 
The  Laplace  transform  is  basically  a  generalized  Fourier  transform  and  the 
z-transform  is  a  generalized  discrete  Fourier  transform  (DFT).  Since  we  will 
eventually  be  using  the  FFT  (an  efficient  computer  algorithm  for  computing  the 
DFT)  to  perform  the  necessary  transformations,  we  will  present  the  following 
discussion  in  terms  of  the  z-transform.  Usually,  the  z-transform  is  defined 
for  sequences  or  sampled  functions  rather  than  for  continuous  functions,  but, 
for  our  purposes,  we  will  extend  the  definition  to  continuous  functions. 


Let  us  define  our  z-transform  as  follows  (compare  to  equation  (A-4)  in 
appendix  A) , 


which  reduces  to  the  Fourier  transform  if  the  complex  variable  z  is  replaced 
by  eit.  Notice  that  the  transform  is  written  in  terms  of  w  instead  of  z  since 
we  will  forward  transform  the  log-spectrum  to  enter  the  cepstrum  domain. 

Figure  4  shows  the  z-plane  and  the  circular  contour  (a  unit-radius  circle 
centered  on  the  origin)  followed  in  the  special  case  of  the  Fourier  transform. 
As  real  to  increases,  the  integration  progresses  around  this  circle  in  the 
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Figure  4.  Fourier  transform  path  in  complex  z-plane. 
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counter-clockwise  direction.  In  order  to  evaluate  the  Fourier  transform,  we 
only  need  to  know  the  values  of  the  spectrum  along  this  circular  path.  It  is, 
however,  useful  to  examine  the  spectrum  in  the  z-plane  near  this  path  as  its 
behavior  near  the  path  affects  its  value  on  the  path.  For  example,  we 
recognize  the  difficulty  involved  when  a  zero  of  the  spectrum  falls  on  the 
Fourier  transform  path:  the  complex  logarithm  becomes  indeterminate.  If  a 
zero  is  near  the  path,  the  trouble  is  not  as  obvious  but  a  spectral  minimum 
will  result  accompanied  by  a  rapid  phase  change. 

To  examine  the  influence  of  zeros  close  to  the  unit  circle  in  the  z-plane, 
let  us  consider  a  spectrum  defined  by  a  single  zero  and  no  poles, 

Fo(z)  =  Ao(z-zo) 

where  F0(elw)  =  F  (z)  for  z  on  the  unit  circle.  This  situation  is  pictured 
in  figure  5  along  with  two  vectors  corresponding  to  the  difference  (z-zQ) 
for  two  values  of  z  on  the  unit  circle.  As  the  frequency  increases  from 
ul  to  w2,  the  length  of  the  vector  from  z  to  z  decreases  to  a  minimum  at  the 
point  of  closest-approach  of  the  circle  t8  z0  and  then  increases  again.  For 
the  same  progression  of  frequency,  the  phase  continually  decreases  (clockwise 
vector  rotation)  slowly  at  first,  then  rapidly  past  the  closest-approach  point, 
then  slowly  again.  These  changes  are  illustrated  in  figure  6.  If  the  zero 
is  inside  the  unit  circle,  the  variations  are  similar  but  the  phase  increases 
with  increasing  frequency. 

Poles  will  appear  in  the  z-plane  if  there  is  any  trapped  propagation  in 
the  ocean  floor  (that  is,  ducting  by  sediment  layering  or  gradients)  or  one  of 

several  types  of  boundary  waves.  These  effects  are  discussed  in  detail  in  the 

companion  volumel  to  this  report;  but,  in  short,  poles  affect  the  phase  in  a 

manner  similar  to  zeros.  When  a  pole  is  near  the  unit  circle  the  phase  will 

also  change  rapidly  as  the  contour  passes  by  the  pole.  In  the  case  of  poles, 
we  also  know  from  stability  theory  that  for  physical  (naturally  stable)  system, 
the  poles  will  be  inside  the  unit  circle  and  so  we  can  predict  what  direction 
the  induced  phase  change  will  be.  Since  the  factor  (z  -  z0)  is  in  the  denomina¬ 
tor  for  a  simple  pole,  the  direction  of  the  phase  change  is  opposite  to  that 
associated  with  a  zero.  Hence,  the  phase  will  decrease  rapidly  near  a  pole 
just  inside  the  unit  circle. 

Since,  in  practice,  the  spectrum  is  sampled,  there  often  is  not  enough 
information  to  establish  the  direction  of  these  very  rapid  phase  changes. 

For  a  zero  very  close  to  the  unit  circle,  the  phase  change  is  roughly  ±tt  and 
this  is  the  same  magnitude  as  the  ambiguity  in  phase  of  a  complex  number. 

(In  other  words,  if  we  guess  that  the  change  was  tt  when  it  actually  was  -it, 
the  error  is  2tt  which  cannot  be  resolved  from  the  complex  number.)  As  we  will 
see  in  the  next  section,  these  phase  changes  can  be  so  rapid  for  real  signals 
that  even  by  increasing  the  sampling  rate  of  the  spectrum  they  are  frequently 
misinterpreted.  In  short,  the  sampled  spectrum  itself  is  adequately  determined 
but  the  presence  of  zeros  (or  poles)  in  the  spectrum  and  the  requirement  for 
absolute  phase  make  the  log-spectrum  very  difficult  to  compute. 
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Figure  5.  Construction  of  response  for  a  single  zero  close  to  the  Fourier  transform  path. 
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Figure  6.  Magnitude  and  phase  changes  in  spectrum  for  transform  path  that 
passes  near  a  zero  in  the  z-plane. 
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In  general,  the  spectrum  comprises  a  large  number  of  zeros  in  the  neighbor¬ 
hood  of  the  unit  circle  and  it  is  sometimes  useful  to  move  these  zeros.  For 
example,  if  there  were  several  zeros  very  close  to  the  unit  circle,  we  could 
simplify  the  phase  reconstruction  if  we  could  move  them  away  from  the  circle. 

This  is  accomplished  by  exponentially  weighting  the  original  time  series, ' 
as  follows, 

g(t)  =  f(t)at.  (9) 

The  resultant  spectrum  is, 

.  /»  oo 

G  ( e1  Ll>)  =  J  f(t)aV1a,tdt  (10) 

0 

=  F(ei(7a). 

For  a  less  than  one,  this  is  equivalent  to  sampling  the  spectrum  along  a  circle 
outside  of  but  concentric  with  the  unit  circle.  From  another  point  of  view, 
this  operation  is  equivalent  to  shrinking  the  spectrum  symmetrically  around  the 
z-plane  origin  and  still  transforming  along  the  unit  circle.  By  using  this 
exponential  weighting,  we  can  shift  zeros  that  are  close  to  the  unit  circle 
inside  and  farther  away  from  the  circle. 

This  operation  of  exponential  weighting  can  help  in  the  reconstruction  of 
phase  but  it  must  be  applied  sparingly.  The  form  of  the  factor  at  that  multi¬ 
plies  the  time  series  is  such  that,  even  for  a  close  to  one,  the  latter  portions 
of  the  time  series  can  be  suppressed  substantially.  Fortunately,  very  small 
amounts  of  exponential  weighting  (0.995  <  a  <  1.00)  can  significantly  improve 
the  reconstruction  process. 5 

To  conclude  this  discussion  of  deconvolution  theory,  we  should  mention 
the  concept  of  minimum-phase. 8  In  general,  both  the  real  and  the  imaginary 
parts  of  the  spectrum  are  needed  to  reconstruct  the  original  time  series,  but, 
if  the  time  series  is  causal  (that  is,  zero  for  time  less  than  zero),  then  the 
original  series  and,  therefore,  the  imaginary  part  of  the  spectrum  can  be 
recreated  from  the  real  part  of  the  spectrum.  Thus,  the  real  and  imaginary 
parts  of  the  spectrum  of  a  causal  time  signal  are  related.  In  particular,  the 
imaginary  part  is  the  Hilbert  transform  of  the  real  part.  This  transform  is 
very  simple  to  apply:  the  inverse  Fourier  transform  of  the  spectrum  with  only 
the  real  part  nonzero  is  taken,  the  resulting  time  series  is  set  to  zero  for 
all  time  less  than  zero,  and  then  the  forward  Fourier  transform  is  taken. 

The  result  is  a  spectrum  with  the  real  part  unchanged 'and  the  appropriate 

imaginary  part. 

This  same  concept  can  be  extended  to  the  log-spectrum.  In  this  case,  the 
real  part  (the  log-magnitude  of  the  spectrum)  is  well-defined  and  single-valued. 

It  would  be  convenient  if  we  could  extract  the  imaginary  part  (the  phase  of 
the  spectrum)  from  the  real  part.  This  is  possible  if  the  cepstrum  is  zero 

for  time  less  than  zero.  This  is  the  condition  for  minimum-phase  and,  if  it 

is  true  for  a  signal,  the  imaginary  part  of  the  log-spectrum  can  be  calculated 
from  the  Hilbert  transform  of  the  log-magnitude  spectrum.  Under  certain 
circumstances,  we  will  be  able  to  use  this  technique  to  eliminate  the  phase 
reconstruction  problem;  however,  the  minimum-phase  condition  is  quite  restrictive. 
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Oppenheim  and  Schafer8  show  that  if  all  of  the  poles  and  zeros  of  a 
spectrum  are  inside  the  unit  circle,  then  the  signal  is  minimum-phase.  (Only 
the  poles  are  required  to  be  inside  the  circle  for  a  stable  signal.)  As  we 
have  seen,  we  can  force  all  of  the  zeros  inside  the  unit  circle  by  sufficient 
exponential  weighting;  however,  this  can  lead  to  obliteration  of  most  of  the 
time  series.  Hence,  we  cannot  indiscriminately  force  a  signal  to  be  minimum- 
phase  and  expect  to  get  meaningful  results.  If  a  single  convolutional  component 
of  the  signal  is  almost  minimum-phase  and  the  other  components  can  be  averaged 
out  by  summing  many  log-spectra,  we  can  take  advantage  of  a  small  amount  of  - 
exponential  weighting  and  the  Hilbert  transform  to  extract  the  nearly  minimum- 
phase  component.  We  will  outline  an  application  for  this  procedure  in  the 
next  section. 

Before  leaving  the  subject  of  minimum-phase,  consider  the  type  of  signal 
that  is  minimum  phase.  Since  we  know  that  exponential  weighting,  in  effect, 
moves  zeros  (and  poles)  toward  the  origin  (for  a  <  1),  and  we  know  that  a 
signal  is  minimum-phase  if  all  of  its  poles  and  zeros  are  within  the  unit 
circle,  we  would  expect  that  a  time  series  that,  on  the  average,  decays  ex¬ 
ponentially  would  be  minimum-phase.  This  is  generally  true.  The  waveform  of 
an  explosion  with  its  high-level  shock  front  and  decreasing  amplitude  bubble 
pulses  is  nearly  minimum-phase.  A  series  of  reflections  of  decreasing  amplitude 
from  a  normal  incidence  reflection  from  a  layered  medium  would  be  almost  min¬ 
imum-phase.  On  the  other  hand,  a  near-grazing  reflection  from  the  ocean  floor 
with  refracted  arrivals  appearinq  later  and  stronger  than  directly  reflected 
arrivals  would  not  be  minimum-phase.  This  intuitive  view  of  minimum-phase  is 
of  some  value  in  deciding  on  the  method  for  phase  reconstruction,  but  it  cer¬ 
tainly  should  not  be  accepted  as  proof  of  the  nature  of  any  signal.  Since 
the  process  described  in  this  report  is  nonlinear,  each  type  of  signal  will 
have  to  be  individually  examined  for  peculiarities.  In  the  nonlinear  world, 
generalizations  are  dangerous. 


APPLIED  HOMOMORPHIC  DECONVOLUTION 


In  this  section,  we  will  consider  the  problems  and  techniques  associated 
with  the  application  of  deconvolution  theory  to  real  signals.  These  techniques 
tend  to  be  application-dependent  and  so  some  modifications  will  undoubtedly 
be  necessary  in  order  for  these  processes  to  work  on  types  of  signals  other  than 
those  considered  herein.  We  have  concentrated  on  deconvolution  of  seismic 
measurements  of  the  ocean  floor  made  with  explosives  and  sonobuoys  over  wide 
ranges  of  incident  angle.  As  a  result,  the  following  discussion  will  be  tailored 
toward  these  types  of  signals;  however,  the  principles  involved  should  help 
the  designer  of  homomorphic  processors  for  other  types  of  convolved  signals. 

Before  we  describe  the  mechanics  of  a  homomorphic  deconvolution  process, 
let  us  review  the  major  obstacles.  For  our  purposes,  the  desired  end  product 
of  deconvolution  is  the  system  impulse  response  (that  is  the  received  signal  for 
a  sharp  excitation  pulse  after  reflection  from  the  ocean  floor).  This  response 
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is  the  system  output  with  all  frequencies  in  the  band  of  interest  excited. 
Unfortunately,  the  explosion  does  not  excite  all  frequencies  uniformly. 
Furthermore,  filtering  in  the  data  acquisition  system  eliminates  still  more 
of  the  desired  information.  The  source  waveshape  spectrum  rolls  off  quickly 
below  the  fundamental  frequency  of  the  bubble  oscillation  (approximately  50  Hz 
for  a  2.2  lb  TNT  charge  at  800  ft)  and  above  that  there  are  many  nulls  in  the 
spectrum.  In  spite  of  the  very  high  signal -to-noise  ratio  signal  produced  by 
the  explosion,  the  noise  in  these  nulls  can  dominate.  This  information  is 
lost  in  any  kind  of  deconvolution  process. 

Inadequate  sampling  of  the  spectrum  creates  another  problem  unique  to 
homomorphic  deconvolution.  While  the  signal  is  adequately  sampled  for  linear 
operations,  the  phase  spectrum  can  be  grossly  undersampled.  Consider  the 
sampled  spectrum  of  an  actual  measurement  (with  noise)  of  which  four  adjacent 
spectrum  points  are  plotted  in  figure  7.  This  is  a  properly  sampled  signal 
according  to  Nyquist's  criterion  but  there  is  no  clear  connection  from  one 
point  to  the  next.  Although  it  becomes  expensive  quickly,  we  can  increase 
the  sampling  of  the  spectrum  by  adding  zeros  to  the  end  of  the  time  series 
prior  to  transforming.  For  example,  if  the  original  time  series  is  256  points 
long  and  256  zero  points  are  appended,  the  resulting  spectrum  will  still  cover 
the  same  frequency  range  but  will  be  sampled  twice  as  densely. 

Figure  8  shows  the  same  four  points  as  figure  7  but  with  a  curve  super¬ 
imposed  by  increasing  the  spectral  sampling  by  a  factor  of  16.  Now  the 
sequence  of  phase  (and  magnitude)  from  point  to  point  is  clear  and  quite  un¬ 
expected.  The  phase  change  from  one  point  to  the  next  is  very  large  (the  order 
of  2 tt )  and  at  one  point  the  curve  passes  so  close  to  the  origin  that  this 
large  increase  in  resolution  was  necessary  for  establishing  the  path.  If  we 
pass  the  origin  on  the  wrong  side,  we  permanently  accumulate  a  2ir  error  in  the 
phase  of  all  succeeding  points.  This  close  pass  is  the  result  of  a  zero  very 
near  the  unit  circle  in  the  z-plane.  By  the  direction  of  the  phase  change 
(a  decrease),  we  know  that  this  zero  must  be  outside  the  unit  circle  and  we 
cannot,  therefore,  be  dealing  with  a  minimum-phase  process. 

Some  attempts  have  been  made  to  reconstruct  phase  by  integrating  the  phase 
derivative. 5  This  sounds  round-about,  but  the  derivative  of  phase  with  respect 
to  frequency  is  single-valued,  whereas  the  phase  itself  is  not.  The  phase  is, 


<(>  =  tan"1  (y/x) 


(11) 


where 


G  =  x  +  iy,  the  complex  spectrum  point, 


and  this  can  only  be  computed  as  a  principal  value  without  additional  informa¬ 
tion  about  the  absolute  phase.  On  the  other  hand,  the  derivative. 


d<j> 

do) 


i  ± 

1  +  (y/x)2  dw 


=  xy'  -yx 
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Figure  7.  Four  adjacent  spectrum  points  of  a  properly-sampled  real  signal 
plotted  in  the  complex  spectrum  plane. 
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Figure  8.  Actual  spectrum  progression  between  the  four  points  of  figure  7. 
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where  the  primes  indicate  differentiation  with  respect  to  to,  is  single-valued. 
Given  a  suitable  starting  value,  this  expression  can  be  integrated  to  give 
absolute  phase.  Unfortunately,  this  method  is  extremely  susceptible  to  under¬ 
sampling.  The  phase  derivative  curve  and  our  four  sample  points  are  shown  in 
figure  9.  Some  of  the  peaks  in  the  derivative  curve  are  missed  completely  and 
even  when  one  is  hit  the  behavior  of  the  curve  is  very  poorly  represented. 

Thus,  the  integration  produces  erroneous  values.  Tribolet^  attempted  to  solve 
this  problem  with  an  adaptive  scheme  that  compares  the  integrated  phase's 
principal  value  with  the  value  computed  by  the  arc-tangent.  This  procedure  .is 
somewhat  more  reliable,  but  when  a  sample  point  lands  on  one  of  the  phase 
derivative  peaks  the  phase  can  be  off  by  one  or  more  complete  cycles.  This 
happened  frequently  enough  to  abandon  the  phase  derivative  method. 

Since  the  spectrum  is  undersampled,  any  reconstruction  method  must  first 
increase  the  sampling  density  in  the  frequency  domain.  Large  increases  in 
resolution  by  zero-filling  the  time  series  are  not  economical,  but  this  is  a 
good  technique  to  use  in  moderation.  The  smooth  trajectory  of  the  curve  in 
figure  8  suggests  that  a  curve-fit  might  be  successfully  applied  to  a  moderate 
sampling  increase.  First,  the  spectrum  sampling  rate  is  increased  by  a  factor 
of  eight.  Then  these  points  are  fit  with  a  series  of  complex  cubic  splines 
based  on  sets  of  six  adjacent  points.  Several  curve-fits  were  evaluated  and 
this  seemed  to  work  the  best  (see  figure  10);  however,  even  the  six-point  spline 
sometimes  passed  on  the  wrong  side  of  the  origin. 

At  this  point,  a  decision  must  be  made  as  to  how  to  proceed.  If  we  are 
content  that  each  point  of  the  spectrum  be  considered  valid  (that  is,  not 
corrupted  by  noise),  we  can  develop  an  algorithm  (described  below)  to  locate 
the  point  of  closest-approach  of  the  spectral  trajectory  to  the  origin.  Having 
found  the  frequency  corresponding  to  this  point,  we  can  calculate  the  DFT  at 
that  point  and  use  that  point  to  resolve  the  path  ambiguity.  Since  there  are 
only  a  limited  number  of  very  close  passes  in  a  typical  signal,  the  discrete 
transform  calculations  are  not  expensive.  This  approach  does  presume  that  the 
spectrum  value  is  a  valid  signal  point  even  though  it  is  of  very  low  level. 

We  would  expect,  at  least  in  some  cases,  that  noise  would  dominate  in  these 
spectral  nulls  and  we  do  not  want  to  permit  noise  to  drive  the  reconstruction 
process. 

A  procedure  which  has  not  been  tested  adequately  is  based  on  a  weighted 
curve-fit  of  the  spectrum  points.  The  spectrum  is  computed  for  increased 
resolution  (four  times  is  probably  sufficient)  and,  as  before,  consecutive  sets 
of  the  points  are  fit  but  in  a  least-square  sense  with  the  higher  amplitude 
values  weighted  more  than  the  lower  values.  Some  type  of  bridging  scheme  like 
this  has  potential  for  enhancing  phase  reconstruction  in  the  presence  of  noise. 

The  physics  of  the  signal  and  the  ocean  medium  may  also  be  used  to  provide 
some  constraints  on  the  phase  reconstruction  process.  Some  preliminary  studies 
of  the  phase  progression  in  real  systems  have  been  summarized  in  the  companion 
volumel  to  this  report,  and  we  have  already  discussed  some  of  these  results  in 
connection  with  poles  in  the  spectrum.  We  can,  at  least,  tell  beforehand  which 
direction  the  phase  change  will  be  in  near  system  resonances  connected  with 
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Figure  9. 


Sampling  of  phase  and  phase  derivative  for  a  typical  real 
signal  sampled  at  the  Nyquist  rate. 
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trapped  or  ducted  propagation  and  boundary  waves.  Considerably  more  work 
needs  to  be  done  in  this  area  in  order  to  better  understand  the  behavior  of 
the  phase  spectrum  of  transient  signals. 

The  technique  used  for  most  of  the  work  presented  here  is  the  procedure 
based  on  locating  the  points  of  closest-approach  and  computing  the  DFT  at  this 
point.  Several  steps  are  involved  here:  first,  close  pass  situations  must  be 
detected;  next,  the  discrete  transform  must  be  evaluated;  finally,  the  path 
route  around  the  origin  must  be  chosen  based  on  this  calculation. 

In  the  interests  of  computational  speed,  a  fairly  crude  (and  conservative) 
scheme  is  used  to  detect  potential  close-pass  situations.  If  the  angular 
difference  between  successive  points  of  the  "oversampled"  spectrum  is  less  than 
90  ,  we  assume  that  the  spectral  trajectory  between  these  points  can  be  ade¬ 
quately  represented  by  a  straight  line.  By  considering  each  of  the  two  points 
as  a  vector  emanating  from  the  origin  and  taking  the  dot  product  of  these 
vectors,  the  sign  of  the  cosine  of  the  included  angle  can  be  evaluated  very 
rapidly  since,  if. 


A  =  (xl5  yi) 
b  -  (x2,  y 2) 

then 


/A//B/cos<j>  =  A  •  B  =  xxx2  +  y.y„  (12) 

where  A,  B  are  vectors  in  the  spectrum-plane.  If  the  dot  product  is  positive, 
then  the  included  angle  <j>  must  be  less  than  90°  and  we  can  proceed  directly  to 
detection  of  path  direction. 

On  the  other  hand,  if  the  included  angle  is  greater  than  90°,  we  fit  the 
set  of  six  adjacent  points  centered  on  the  trouble  spot  with  a  complex  spline 
(actually  two  real  splines,  one  on  the  real  parts  of  the  points  and  one  on  the 
imaginary  parts).  A  fast,  bisection  search  is  then  performed  to  locate  the 
closest-approach  point  and  the  DFT  is  computed  at  the  corresponding  frequency. 
The  DFT  need  not,  however,  be  taken  on  the  entire  time  series,  which  is  mostly 
zero  (having  been  padded  with  zeros  to  increase  resolution).  Instead,  we  take 
the  DFT  of  the  unpadded  time  series  and  then  adjust  the  phase  of  the  spectrum 
point  to  reflect  the  shift  of  the  original  series  within  the  field  of  zeros. 
The  path  direction  is  then  determined  from  this  intermediate  point. 

One  important  feature  of  the  phase  tracking  procedure  is  that  the  actual 
path  of  the  spectrum  is  immaterial  as  long  as  a  count  is  maintained  of  the  net 
number  of  orbits  of  the  spectrum  origin.  With  this  count  and  the  principal 
value  of  the  phase  at  a  point,  the  point's  absolute  phase  can  be  determined 
merely  by  adding  or  subtracting  a  corresponding  multiple  of  2 u.  The  count  of 
origin  orbits  can  be  kept  by  counting  the  net  number  of  crossings  of  the  nega¬ 
tive  real  axis  in  the  counter-clockwise  direction.  This  implies  a  branch  at 
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±tt  which  is  consistent  with  the  principal  value  branch  selected  by  the  com¬ 
puter  arc-tangent  function.  When  it  is  necessary  to  compute  a  close-approach 
point,  the  crossing  detection  is  done  in  two  stages:  from  one  spectrum  point 
to  the  approach  point  and  then  from  the  approach  point  to  the  next  spectrum 
point. 

We  have,  in  essence,  solved  the  phase  reconstruction  problem  for  a  sampled 
signal.  In  this  implementation  it  is  still  possible  for  errors  to  occur,  so 
as  interative  test  procedure  is  used.  Since  we  must  remove  the  linear  (ramp) 
component  of  the  phase  change  from  one  end  of  the  spectrum  to  the  other  by 
shifting  the  time  series,  we  can  redo  the  phase  reconstruction  on  this  shifted 
series.  The  results  should  be  the  same  after  correcting  for  the  linear  phase 
shift  induced  by  the  time  shift.  The  relationship  between  these  two  shifts 
can  be  found  by  taking  the  Fourier  transform  of  a  shifted  time  series, 

03 

/f(t-t0)e-i«tdt  =  eiajt0  (13) 

o 

where 


F(uj)  =  Fourier  transform  of  f(t). 

The  first  spectrum  point  is  the  zero-frequency  point  so  the  absolute  phase 
there  is  set  to  zero.  The  highest  point  in  the  spectrum  u>t,  after  the  initial 
phase  reconstruction,  has  some  absolute  phase  <f>t  that  can  be  removed  by  time 
shifting  according  to  equation  (13)  where 


or 


tQ  =  -*t/V  (14) 

This  will  also  affect  every  other  point  in  the  spectrum  as  follows, 

Fs<“>  3  (16) 

where  Fs(u>)  =  Fourier  transform  of  the  shifted  time  series.  Thus,  the  phase 
shift  is  in  the  form  of  a  ramp:  zero  at  zero  frequency,  <j>t  at  the  highest 
frequency,  and  linear  with  frequency  in  between. 

If  we  had  complete  confidence  in  the  initial  phase  reconstruction,  we 
could  remove  the  resulting  phase  ramp  by  a  suitable  time  shift  and  proceed  with 
the  deconvolution.  However,  the  process  is  not  always  error-free,  so  the 
reconstruction  is  repeated  on  the  shifted  time  series.  In  addition  to  the 
usual  test  for  reconstruction  near  zeros  or  poles,  each  new  point  is  compared 
with  each  corresponding  point  in  the  previous  pass  to  see  if  the  phase  has 
changed  proportionately  following  equation  (15).  Any  time  a  discrepancy  is 
found,  the  spline  fit  and  closest-approach  calculations  are  performed  to  resolve 
the  questionable  point. 
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When  this  second  phase  reconstruction  has  been  completed,  the  absolute 
phase  of  the  last  point  in  the  spectrum  should  be  zero  (since  the  phase  ramp 
was  removed  after  the  first  pass).  If  not,  the  new  ramp  is  removed  and  the 

reconstruction  is  performed  again.  In  most  cases,  this  is  sufficient.  If 

the  endpoint  phase  is  still  not  zero  the  procedure  can  be  repeated,  but  this 
may  be  an  indication  that  the  signal  cannot  be  adequately  resolved  by  this 
technique. 

If,  on  the  other  hand,  the  second  or  third  reconstruction  is  successful- 
(endpoint  phase  goes  to  zero),  then  the  absolute  phase  values  can  be  stored 
as  the  imaginary  part  of  the  log-spectrum.  The  real  part  is  formed  by  the 
natural  logaritm  of  the  spectrum  magnitude  and  this  completes  the  calculation 
of  the  log-spectrum.  Since  generation  of  the  log-spectrum  is  the  key  process 
in  this  method  of  deconvolution,  a  listing  of  a  computer  code  for  this  operation 
is  given  in  appendix  B. 

Most  of  this  discussion  so  far  has  been  centered  on  the  phase  reconstruc¬ 
tion  problem  because  this  is  the  most  difficult  aspect  of  homomorphic  decon¬ 
volution.  The  other  processes  cannot,  however,  be  done  carelessly.  Proper 
conditioning  of  the  signal  prior  to  phase  reconstruction  is  vital,  and  this 
conditioning  will  be  discussed  in  some  detail  in  the  next  section.  In  short, 
this  preconditioning  involves  proper  sampling,  filtering,  and  exponential 
weighting. 

Once  the  log-spectrum  has  been  calculated,  a  standard  FFT  is  applied  to 
compute  the  complex  cepstrum  (unless  the  deconvolution  filter  has  been  designed 
to  operate  directly  on  the  log-spectrum).  Then  the  portion  of  the  cepstrum, 
corresponding  to  the  undesired  part  of  the  signal,  is  removed.  Simply  setting 
these  points  to  zero  is  often  sufficient.  For  example,  if  we  want  to  eliminate 
the  bubble  oscillations  of  an  explosive  source,  we  zero  a  portion  of  the  cepstrum 
around  the  bubble  oscillation  period  (that  is,  the  time  between  repetitions). 

We  might  also  zero  a  small  band  of  time  corresponding  to  the  repetition  rate 
generated  by  unwanted  multipath  reflections. 

Having  isolated  the  desired  portion  of  the  cepstrum,  we  then  perform  the 
inverse  FFT  to  return  to  the  log-spectrum  and  then  perform  the  complex  expo¬ 
nential  is  computed  as  follows, 


exp(x+iy)  =  exe^ 

This  function  is  completely  unambiguous. 

Finally,  the  inverse  FFT  is  applied  once  again  and  the  resulting  time 
series  is  exponentially  weighted  by  the  inverse  of  whatever  factor  was  used  to 
weight  the  original  time  series.  After  shifting  the  series  back  to  its  orig¬ 
inal  starting  point  in  time,  we  have  completely  recovered  the  deconvolved 
signal.  An  example  of  the  complete  process  is  given  in  figure  11  where  the 
source  waveshape  (A)  is  removed  from  a  signal  (A*B)  with  multiple  arrivals. 

The  received  signal  was  synthetically  generated  by  convolution  of  a  replica  of 
the  explosion's  bubble  pulse  and  a  three-path  reflector  series  (B).  While  the 
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Figure  11.  Example  of  homomorphic  deconvolution  of  a  received  signal 
(A*B)  made  up  of  a  source  excitation  (A)  and  a  multipath 
medium  response  (B). 
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source  waveform  was  known  to  construct  the  received  waveform,  this  information 
was  not  used  in  the  deconvolution.  The  deconvolution  was  performed  by  zeroing 
the  cepstrum  for  all  time  at  or  below  the  bubble  oscillation  period. 

Additive  noise  in  the  time  domain  can  be  tolerated  up  to  a  certain  point 
without  adversely  affecting  the  process.  Beyond  this  point,  however,  the  phase 
reconstruction  fails  because  the  spectrum  nulls  and  the  attendant  phase  jumps 
become  dominated  by  the  characteristics  of  the  noise.  Hopefully,  the  pre¬ 
viously  mentioned  null-bridging  technique  will  reduce  this  problem.  Until  - 
null-bridging  is  perfected,  this  type  of  deconvolution  is  probable  restricted 
to  high  signal-to-noise  signals.  Fortunately,  most  of  the  ocean  floor  measure¬ 
ments  with  which  we  are  concerned  are  high  signal-to-noise,  thanks  to  the  high 
intensity  sound  generated  by  the  explosive  charges. 

In  the  preceding  section,  the  concept  of  minimum-phase  was  introduced. 

For  reasonably  deep,  underwater  explosions,  the  source  waveshape  is  almost 
minimum-phase.  This  property  can  be  exploited  to  extract  the  source  waveform 
itself.  Since  the  measurements  on  which  we  are  operating  are  taken  over  a 
wide  range  of  angles,  the  only  component  of  the  signal  common  to  all  received 
signals  at  a  given  site  is  the  source  waveshape.  If  we  then  add  a  large  number 
of  the  log-spectra  of  these  measurements,  the  components  related  to  the  ocean 
response  will  tend  to  cancel  while  the  source  components  will  build  up.  We 
can  avoid  the  problem  of  phase  reconstruction  by  sufficient  exponential 
weighting  to  insure  that  the  source  waveform  is  minimum-phase  and  then  by 
Hilbert  transforming  to  generate  the  phase  of  the  averaged  log-spectrum. 10 
The  ocean  response  need  not  be  minimum-phase  because  these  components  are 
averaged  out,  each  measurement  having  a  different  geometry  and  arrival  structure. 

After  the  averaged  log-spectrum  has  been  completed  by  Hilbert  transforma¬ 
tion  of  the  averaged  log-magnitude  spectrum,  the  spectrum  and  corresponding 
time  series  are  computed.  Figure  12  shows  an  example  of  the  isolation  of 
the  source  waveform  from  20  seismic  measurements,  each  with  a  complicated 
arrival  structure.  The  source  waveshape  extracted  is  very  similar  to  the  ex¬ 
pected  explosion  waveshape.il 

Once  the  source  waveshape  is  established  in  this  way,  one  of  the  more 
conventional  deconvolution  processes  (such  as  inverse  filtering)  may  be  used 
to  obtain  the  ocean  floor  response.  One  practical  point  worth  mentioning  here 
is  that  the  source  waveform  of  each  measurement  may,  in  fact,  be  slightly 
different  because  the  charge  does  not  always  detonate  at  the  same  depth. 

The  depth  of  the  explosion  (along  with  the  relatively  constant  charge  weight) 
determines  -the  period  of  the  bubble  oscillation.  We  can,  however,  measure 
this  period  quite  accurately  by  the  associated  peak  in  the  cepstrum  and  correct 
each  measurement  to  the  same  depth  by  stretching  or  shrinking  the  time  axis 
prior'  to  averaging.  The  cepstrum  computed  from  the  log-magnitude  with  zero 
phase  is  adequate  for  measuring  the  bubble  period. 

It  is  also  possible  to  perform  this  averaging  deconvolution  to  extract 
the  ocean  floor  response  if  the  geometry  is  such  that  the  reflection  response 
is  almost  minimum-phase  and  the  source  waveshape  can  be  forced  to  vary  from 
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Figure  12.  Reconstruction  of  explosive  source  waveform  by  log-spectral  averaging. 
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measurement  to  measurement.  A  fairly  high  angle  (i.e.,  near-vertical) 
measurement  may  satisfy  the  minimum-phase  requirement  with  slight  exponential 
weighting  and  several  shots  of  varying  depths  and  sizes  could  be  used  to  vary 
the  source  waveshape. 

Also,  it  should  be  possible  to  reduce  the  noise  problem  for  any  type  of 
deconvolution  by  averaging  spectra  or  log-spectra  of  a  number  of  measurements 
with  the  same  geometry  and  source  waveshape.  Neither  of  the  last  two  proposals 
has  been  tried  as  each  would  require  a  new  experimental,  design  (probably  in 
the  form  of  a  shipboard  experiment),  but  they  have  enough  promise  to  encourage 
their  consideration  in  future  experimentation. 


SIGNAL  CONDITIONING 


As  we  have  emphasized,  the  success  of  this  signal  process  is  dependent  on 
the  quality  of  the  signal  on  which  the  process  operates.  Such  general  and 
powerful  principles  as  superposition  or  Nyquist's  sampling  theorem  have  no 
equivalents  in  nonlinear  processing.  Each  type  of  signal  must  be  treated  on 
its  own  merits  and  must  be  conditioned  properly  in  order  to  insure  a  valid 
deconvolution.  In  this  section,  we  will  outline  some  of  the  ways  in  which 
signals  should  be  conditioned  as  a  part  of  homomorphic  deconvolution. 

First,  let  us  review  the  basic  problems  with  nonlinear  deconvolution. 

In  deconvolution,  we  are  trying  to  reconstruct  the  response  of  a  system  at  all 
frequencies  (the  impulse  response)  from  the  system  output  for  an  excitation 
that  does  not  contain  all  frequencies.  This  means  that,  for  those  portions  of 
the  spectrum  for  which  we  have  little  signal,  noise  will  corrupt  the  convolu¬ 
tion.  In  addition,  for  operations  with  the  log-spectrum,  minimum-phase  signals 
are  ideal;  however,  we  usually  must  deal  with  mixed-phase  signals. 

To  avoid  processing  portions  of  the  spectrum  that  are  not  well -excited  by 
the  source,  we  should  filter  and  translate  the  original  received  signal  for 
maximum  energy  across  the  remaining  band.  Remember  that  we  are  not  free  to 
select  any  band  for  the  impulse  response;  the  spectrum  of  the  source  determines 
what  regions  of  the  frequency  response  can  be  found.  In  general,  we  will  have 
to  low-pass  filter  to  remove  regions  of  high  frequency  noise  and  then  translate 
to  remove  the  poorly-excited  frequency  range  below  the  explosion's  bubble 
oscillation  frequency.  Some  work  has  been  done  in  breaking  up  the  spectrum 
between  these  1 i m i t s 1 2  so  as  to  eliminate  intervening  spectral  nulls  resulting 
from  the  source,  but  we  will  not  consider  these  techniques.  For  the  present, 
it  is  sufficient  to  drop  the  high  and  low  sides  of  the  spectrum  where  the 
source  energy  is  rolling  off.  The  objective  of  this  first  step  is  to  minimize 
the  number  and  extent  of  low-level  areas  in  the  spectrum. 

Next,  the  signal  should  be  sampled  as  close  to  the  Nyquist  rate  (that  is, 
twice  the  frequency  of  the  highest  significant  component  of  the  signal)  as 
possible.  Undersampling,  of  course,  will  introduce  aliasing.  Oversampling 
may  be  just  as  bad,  though,  as  oversampling  effectively  pushes  the  high 
frequency  limit  of  the  analysis  up  beyond  the  roll-off  of  significant  energy 
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in  the  spectrum.  This  generates  a  region  of  low  spectrum  level  that  will 
be  difficult  to  resolve  during  the  phase  reconstruction.  Consequently, 
filtering  of  the  signal  should  be  done  with  a  fairly  sharp-edged  filter  so 
that  aliasing  can  be  avoided  while  still  maintaining  a  near-optimum  sampling 
rate. 


In  the  type  of  signals  considered  in  this  report,  there  are  generally  a 
large  number  of  zeros  very  close  to  the  unit  circle.  Many  of  these  result 
from  the  nearly-minimum-phase  source  waveform  and  they  are  all  seen  as  large 
dips  in  the  spectrum.  The  phase  reconstruction  can  be  simplified  somewhat 
by  a  small  amount  of  exponential  weighting.  As  we  have  seen,  this  weighting 
moves  the  zeros  toward  the  origin  of  the  z-plane  and,  therefore,  away  from  the 
unit  circle.  Since  the  amplitude  of  the  time  series  decreases  exponentially 
with  time  for  this  weighting,  we  must  first  shift  the  time  series  of  the 
received  signal  toward  the  origin  (t=0)  so  that  the  signal  starts  immediately 
on  the  time  record.  Also,  the  a  in  equation  (9)  should  not  be  less  than 
about  0.99  (and,  of  course,  not  greater  than  1.00).  Below  this  value,  the 
time  series  is  seriously  attenuated  toward  the  end  of  the  record.  As  a  result, 
it  would  be  unusual  if  this  limited  magnitude  of  exponential  weighting  made 
the  entire  signal  minimum-phase,  but  it  may  be  enough  to  make  the  source  com¬ 
ponent  minimum-phase. 

After  the  exponential  weighting  has  been  done,  it  is  often  helpful  to 
shift  the  time  series  again.  Since  the  phase  ramp  in  the  spectrum  is  undesir¬ 
able  and  corresponds  to  a  uniform  shift  or  delay  in  the  time  series,  we  can 
reduce  this  ramp  considerably  by  shifting  the  weighted  time  series  so  that  it 
is  approximately  centered  on  t=0  (equal  energy  on  both  sides  of  t=0,  for 
example).  This  will  not,  of  course,  completely  eliminate  the  phase  ramp,  but 
it  will  simplify  the  first  pass  of  the  phase  reconstruction. 

In  general,  it  is  helpful  to  know  as  much  as  possible  about  the  physics 
involved  in  a  particular  signal  formation  process.  The  physics  of  a  process 
can  indicate  what  convolutional  components  are  likely  to  be  present.  In  ocean 
floor  seismic  measurements,  these  components  may  include  multipath  effects, 
attenuation,  and  phase  shifts.  The  medium-induced  phase  shifts  may  be  particu¬ 
larly  important  as  they  can  introduce  signal  distortion  at  ray  turning  points 
or  at  reflections  from  impedance  discontinuities.  Knowledge  of  the  explosion's 
bubble-pulse  formation  and  oscillation  allows  us  to  guess  that  this  is  close 
to  a  minimum-phase  process  and  permits  us  to  exploit  some  of  the  special  prop¬ 
erties  of  these  processes. 


CONCLUSIONS 


One  of  the  most  important  results  of  this  study,  aside  from  the  develop¬ 
ment  of  the  deconvolution  process,  is  that  the  full  homomorphic  deconvolution 
with  phase  reconstruction  is  highly  dependent  on  the  quality  and  conditioning 
of  the  original  signal.  The  process  works  well  on  signals  that  have  a  high 
s i gnal - to-noi se  ratio,  are  strongly  excited  over  most  of  the  frequency  range 
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of  interest,  and  are  properly  sampled  and  weighted.  While  the  phase  recon¬ 
struction  method  described  above  is  effective  in  many  cases,  sufficient 
degradation  in  the  signal  quality  will  eventually  produce  meaningless  phase 
spectra.  Consequently,  applications  of  this  deconvolution  process  will  require 
careful  checking  of  the  quality  of  the  input  signals. 

This  is  not  to  say  that  the  process  must  be  abandoned  in  the  case  of  less 
strongly  excited  systems;  however,  further  refinements  may  become  necessary. 

The  null-bridging  technique  may  provide  a  means  for  resolving  the  phase  over 
small  regions  of  low  signal  level.  In  addition,  physical  constraints  such  as 
the  phase  change  near  poles  of  the  spectrum  or  at  ray  turning  points  will  be 
helpful  in  predetermining  the  phase  progression.  Development  of  these  tech¬ 
niques  will,  however,  require  substantial  testing  and  experience  with  real 
signals . 

Extraction  of  almost-minimum-phase  components  by  averaging  deconvolution 
is  a  very  effective  process.  In  the  measurements  with  which  we  are  concerned, 
the  only  component  that  this  is  well -suited  to  is  the  source  waveform.  Informa¬ 
tion  about  the  source  waveform  is,  however,  necessary  in  applying  conventional 
deconvolution  and  so  a  hybrid  process  may  be  worthwhile.  Such  a  process  would 
isolate  the  source  waveform  by  log-spectral  averaging  and  then  use  this  source 
waveform  to  develop  an  inverse  filter.  This  would,  of  course,  lose  the  flexi¬ 
bility  of  the  fully  homomorphic  technique  in  that  only  the  impulse  response 
for  the  entire  ocean  system  could  be  found. 

While  tests  on  small  sections  of  experimental  data  have  demonstrated  at 
least  a  limited  value  in  homomorphic  deconvolution,  much  more  experience  must 
be  gained  in  order  to  ascertain  the  extent  of  practical  applications.  This  is 
true  because  the  log-spectrum  and  the  cepstrum  are  not  well-known  domains. 

Until  the  nature  of  signals  in  these  domains  becomes  familar  enough  to  develop 
intuitive  expectations  that  guide  much  of  linear  signal  processing  applications 
day-to-day,  homomorphic  deconvolution  will  remain  a  technique  accessible  to 
only  a  few  specialists.  Hopefully,  much  of  this  experience  will  be  gained  as 
this  process  is  applied  to  the  large  volume  of  seismic  measurements  resident 
at  NAVAIRDEVCEN. 
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TRANSFORM  RELATIONSHIPS 


The  integral  transform  pairs  with  which  we  are  concerned  in  this  report 
are  the  Fourier  transform  pair. 


/oo 

f(t)e"lajtdt 


(A-l) 


f(t) 


—  f 

^  J 


G(w)e1  aitdoj 


and  the  Laplace  transform  pair, 
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which  is  essentially  a  generalized  Fourier  transform.  For  the  special  case  of 
s  =  iu),  the  Laplace  transform  reduces  to  the  Fourier  transform  as  long  as  f(t) 
is  casual  (f(t)  =  o  for  t  <  0)  and  stable  (which  permits  a  =  0). 


Each  of  these  integral  transforms  has  a  discrete  equivalent  that  would  be 
used  on  sampled  series.  The  discrete  equivalent  of  the  Fourier  transform  is 
the  discrete  Fourier  transform  (DFT)  given  by  the  pair, 
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The  discrete  equivalent  to  the  Laplace  transform  is  written  in  terms  of  a  new 
variable  z  which  is  equal  to  es.  This  transform  is  called  the  z-transform  and 
the  corresponding  transform  pair  is, 

oo 

X(z)  =  £  x(n)z”n 

n="~  (A-4) 

x(n)  =  X(z)zn_1dz 

C 

Notice  that  this  transform  is  discrete  in  time  (t  =  nAt)  but  not  in  frequency; 
z  is  a  continuous  variable.  The  contour  C  of  the  inverse  transform  is  any 
closed  contour  lying  entirely  in  the  region  of  convergence  of  the  function  X(z). 
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The  relationship  between  the  various  transforms  can  be  seen  by  comparing 
figures  A-l  and  A-2.  Figure  A-l  shows  the  paths  for  the  Fourier  (F)  and  Laplace 
(L)  inversion  integrals  in  the  z-plance.  The  z-transform  inversion  follows  the 
same  path  as  the  Laplace  inversion.  In  the  case  of  the  DFT,  there  would  be  N 
samples  evenly  spaced  around  the  unit  circle  (the  F  path). 

In  the  s-place  (figure  A-2),  the  circular  contours  unwrap  into  the  verti¬ 
cal  lines  F  and  L  and,  at  the  same  time,  the  entire  z-plane  is  replicated  in 
the  s-plane  as  an  infinite  series  of  horizontal  layers.  The  interior  of  the 
unit  circle  maps  into  many  horizontal  bands  to  the  left  of  the  imaginary  axis 
and  the  exterior  of  the  circle  maps  into  horizontal  bands  to  the  right.  Classi¬ 
cal  stability  theory  says  that  there  may  be  no  poles  in  the  right  half  of  the 
s-plane;  this  corresponds  to  having  all  the  poles  inside  the  unit  circle  in  the 
z-plane.  For  Fourier  analysis,  frequency  starts  at  zero  at  the  origin  of  the 
s-plane  and  increases  vertically  upwards.  In  the  z-plane,  spectral  frequency 
is  zero  at  z  =  +  1  and  increases  counter-clockwise  around  the  unit  circle. 
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Figure  A-l.  Configuration  of  Fourier  and  Laplace  inversion  paths  in  complex  z-plane. 
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s  —  plane 


Figure  A-2.  Configuration  of  Fourier  and  Laplace  inversion  paths  in  complex  s-plane. 
Shaded  area  and  numbered  points  correspond  to  figure  A-l. 
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APPENDIX  B 

COMPUTER  CODE 


B-l 


PROGRAM  TRAC5 ( INPUT , OUTPUT , TAPE20. TAPE40) 
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TRANSFER  256  TIME  POINTS  INTO  LARGER  ARRAY  (FOR 
ZERO-PADDING)  SHIFTED  BY  M  POINTS. 


ADD  OVERFLOW  ON  BOTH  ENDS  OF  DATA  (SEE  STRUCTURE 
IN  COMMON)  TO  ALLOW  SPLINE  FIT  NEAR  ENDS. 


CD 

< 

CL 


in 

in 


O 


cn 

O 


CN 

00 


REPORT  NO.  NADC-82253-30 


00 

CN 

+ 

to 

rr 


in 
0 
< 
q : 


oo  O 
tt  in 

O  O  ^  ^ 


O 


in 

z> 

o 


CL 

CL 


o 

X 


in 
X 
O  I 


00 

2 

o 


G 

CL 


O 

0 


^  CN 

»— •  w 
w  < 

CD  K 
<  <  ~ 
H-  Q  tt 
G  *  • 

o 

CL  —  CN 
CL  I  LL 
LU  h  ^ 


IS) 

< 

X 

a 


CN  CN 

CN 

>S 

CL 

in 

O 

2  CN 

in 

z 

0 

w * 

'w’ 

'w 

X 

1— « 

CL 

0 

3 

<  ^ 

< 

<j 

<  < 

< 

< 

H- 

* 

in 

2 

1— 

•  < 

/— s. 

1  ( 

2 

K  a 

a 

1— 

*-* 

CN 

2 

Ul 

< 

0  H 

2 

a 

0 

<  < 

<S 

< 

3 

w 

H 

O 

L L 

LU 

2  < 

0 

O  O 

Q 

o 

LU 

0 

CL 

CL 

CL 

3 

1 -f 

in 

s. 

<  O 

2 

a 

n  it 

ii 

UJ 

O 

LL 

Ul 

Ul 

X 

LU 

0 

a  * 

UJ 

z 

0 

+ 

II 

CL 

h- 

a 

3 

in 

2 

in 

X 

2 

DC 

< 

CN 

< 

2 

as 

in 

a 

UJ 

< 

t- 

<r 

O  CN 

2 

,-v 

CL 

(-H 

•* 

G 

UJ 

O 

LL 

CN 

i 

<£ 

CN  TT 

CN 

n 

o 

* 

<— s 

-J  2 

b~ 

i 

z 

t— i 

CN 

0 

0 

>w  w 

* 

s—-' 

O 

CL 

CN 

UJ 

M 

< 

1-4 

LU 

Q 

■>“  * 

2 

<1 

<*“ 

CN 

CN 

a 

0 

> 

in 

V 

Ul 

l  CN 

< 

< 

* 

Q 

Q  Q 

a 

a 

a 

O 

0 

< 

0 

CN 

*— i 

■»— 

3 

t— i 

0 

in 

CN 

CN 

<  < 

< 

< 

o 

ll 

< 

CL 

H- 

2 

Ul  O 

X 

1 

O 

H- 

a 

*  < 

< 

G_  CL 

CL 

a 

2 

00 

CL 

CL 

< 

< 

G 

-J 

h— 

H 

t— < 

o 

UJ 

o 

CN  h- 

Q 

2 

2 

< 

< 

X 

h- 

LU 

o 

UJ 

0 

* 

TT 

m 

3 

w  < 

i—i 

0 

| 

* 

00 

o 

0 

CJ 

UJ 

CL 

2 

0 

LU 

CN 

h- 

0  Q 

< 

< 

,  , 

CN 

t*  0  0  0 

UJ 

CN 

LU 

2 

II 

in 

in 

< 

2 

CL 

w 

o 

CL 

0  - 

, — , 

in 

o 

Q 

* 

00 

a 

s 

0. 

< 

< 

O 

*— t 

< 

0 

h» 

o 

2 

a  ^ 

2 

< 

CN 

CJ 

< 

o 

*— 1 

to 

> 

X 

+  o 

LU 

G 

a 

UJ 

CL  ■ 

0 

in  ro 

0 

»— * 

0 

'~v 

X 

lu 

CL 

+ 

CL  « 

in 

CL 

2 

a 

o 

CL  CL 

HH 

•  t 

a 

0 

CO 

1 

< 

a 

f-  <7 

CL 

CD 

o 

CL 

O 

h- 

< 

O 

a 

in  0 

0 

UJ  o 

in 

0  K-4 

0 

< 

a 

in 

3 

UJ 

ll 

LU  Z 

X  o 

X 

u 

0 

2 

CL 

2  * 

DC 

CL 

< 

O  O 

• — ~ 

UJ 

2 

H 

£L 

ii 

*-i 

1—  «-« 

►  H- 

CL 

0 

i 

< 

/— s 

2  CL 

+ 

<  CN 

2 

0 

* 

Q 

i 

CN  CN 

CO 

a 

*— 1 

* 

X 

G  O 

> 

< 

LL 

T— 

<  UJ 

in  w 

0 

CN 

' — ' 

G 

3 

as 

h- 

CD 

in 

in 

an 

-1  CL 

W  O 

UJ 

0 

K-t 

cc 

m  < 

0 

0 

0 

0 

y_^ 

<  < 

< 

< 

* 

— 1 

o 

< 

LL 

CJ 

CN 

* 

■»— 

O 

CN  0 

K 

H- 

in 

2  G  O 

in  lu 

a. 

<  K 

CD 

a 

a 

< 

a 

,  t 

1“  h- 

H- 

K 

< 

O 

2 

*— < 

2 

o 

1 

in  x 

2 

G  u 

LU 

< 

w 

<  X 

UJ 

•  < 

< 

O 

a 

a 

0 

SI- 

<  < 

< 

< 

G 

in 

0 

o 

X 

HH 

CD  H- 

< 

HH 

in 

in 

in 

3  h- 

0 

*-»  Q 

a 

a 

♦— t 

< 

CL 

CN 

o  o 

Q 

Q 

w 

CO 

< 

•  1 

in 

• 

i 

■M 

<  = 

h-  ^ 

O  X 

in 

m 

G 

0 

2 

< 

a  • 

. 

Q 

a 

II  II 

n 

II 

LU 

< 

o 

i-t 

CN 

*— i 

< 

in 

in 

< 

< 

< 

LU  LU 

< 

o  o 

Q 

0 

1 

a 

a 

CJ 

l f) 

2 

2 

in 

It 

W 

w 

UJ 

o 

03 

a 

I- 

CL  h- 

o  o 

0  2 

0 

a 

LU 

< 

LU 

O 

II  W 

M 

< 

a  CL 

II  UJ 

< 

o 

n 

II 

in 

LU  < 

ii 

10  ^ 

# 

2  a. 

< 

0 

— ) 

n 

*  v  ' 

*  v 

*  ' 

X 

a 

*— • 

O 

II 

h- 

G  O 

» 

2 

2 

LU 

X  0 

0  G) 

> _ . 

•a  0 

a 

Lu 

O  u  u 

in 

■*"  CO 

CO 

CL 

X) 

a 

CL  II 

U 

o 

< 

CL  LL 

0  X 

UJ 

< 

LU 

0  1- 

\- 

i-  G 

a 

a 

2  G 

II 

0 

D 

< 

< 

•w  N— * 

a 

< 

CL 

M 

ii 

O  u_ 

Q 

2 

2  G 

* 

Ul 

2 

in 

o 

2 

a 

a  a 

< 

<  O 

ii 

n 

CN 

CN 

—I 

2 

a 

UJ  li- 

X 

in 

t— l 

O  a 

<r  t-H 

< 

in 

3 

< 

Ul 

G 

Ll.  G 

0  2  2 

2 

x  a 

0 

0 

0 

0 

Q  G 

O 

O 

0 

o 

G 

CD  -H 

in 

< 

X 

n 

<J  G 

in  w 

h- 

in 

H- 

LU 

H 

HM  0 

0 

*->  *— t 

CL 

0  m 

a 

0 

_J 

u 

0 

x 

<  <. 

< 

< 

< 

0 

a 

2  X 

CJ 

o 

in 

a 

03  ll 

< 

LU 

M 

G 

u 

0 

2 

DC  DC 

O 

G 

a 

< 

< 

0 

< 

<r 

C3 

Cl  Q. 

CL 

a 

O 

2 

o 

<  in 

2 

o 

< 

X 

2 

< 

CL 

CD 

Q 

O 

< 

< 

a  a 

LU 

in 

0 

0 

2 

0 

0  2 

O 

O 


CL 

cd 

o 

CL 

CL 


0000 


00000 


o 

to 


in 

ID 


O 

r- 


O 

oo 


in 

oo 


O 

a> 


in 

cn 


C 

c 


in 

O 


B-3 


REPORT  NO.  NADC-82253-30 


0 

< 

CL 


o 


in 

in 


r- 

O 


cn 

O 


in 

. 

O 

4-K 

a 

+ 

< 

UJ 

in 

b- 

tn 

UJ 

D 

a 

HH 

Z) 

2  a 

Z 

in 

0 

* 

O 

a 

a 

UJ  < 

< 

•*- 

3 

II 

LU 

O 

CQ 

_ _ , 

a 

u- 

ex 

H> 

in 

\ 

Z 

00  2 

00 

K 

b~ 

0 

CO 

a 

K 

O  0 

in 

in 

• 

in 

3 

a 

< 

o 

o 

o 

U- 

CX  _l 

< 

< 

cn 

Z  • 

0 

a 

2 

2 

a  cq 

CL 

LU 

*-< 

o  o 

Z 

a 

a 

o 

ex 

# 

0  LU 

O 

X 

2 

tn 

n 

a  ex 

in 

< 

> 

a 

< 

o  a 

3 

>  o 

in 

0 

O 

2 

ex 

-J  2 

Z 

< 

tn 

in 

Z  -j 

HH 

LU 

ex 

CX  UJ 

o 

a 

a 

in 

a 

O  — i 

> 

-J 

LU 

UJ  ex 

HH 

a 

a 

a 

x 

3 

< 

LU 

CQ 

CJ 

G. 

a 

< 

Q 

3 

a 

0 

h— 

ex 

O 

z 

o  z 

< 

X 

a 

2 

<  a 

a 

ex 

< 

ex  lu 

in 

2 

< 

X 

1-  o 

CL 

CL  UJ 

z 

a 

> 

0 

X 

3  a 

z 

CJ 

CQ 

UJ 

a 

in 

3 

0 

CJ 

CL 

♦—1 

0 

z 

z 

a 

CO 

a 

a 

o 

a 

Z 

2  0 

0 

< 

UJ  in 

2 

o 

o 

O  3 

a 

< 

tn 

LU 

LU  < 

O 

a 

0 

3 

a 

+ 

CJ  -J 

LU 

HH 

*— 1 

CQ  X 

0 

a 

a 

a 

a 

LU  < 

z 

ex 

CN 

a 

o 

< 

x 

00 

a  > 

*—• 

o 

* 

w 

in  a 

a 

a 

in 

z 

2 

a 

0 

2 

LL 

<  z 

2 

< 

HH 

0 

v 

0 

z 

ex  tn 

ex 

O 

CJ 

X  uj 

< 

* 

< 

in 

a 

O  *-h 

LU 

LU 

LU 

z 

z 

a 

0 

2 

2 

a 

<->  K 

II 

a  I 

b- 

3 

1- 

< 

uj  a 

a 

a 

-  3 

a 

LU 

a 

< 

a 

in  cl 

UJ 

a 

♦“ 

> 

a 

0 

O  0 

2 

CO 

X 

O 

< 

2 

_i 

<  2 

o 

in 

in 

CO  3 

a 

in  2 

a 

0 

lu  on 

> 

nn 

o 

X  O 

o 

O 

< 

0  0 

a 

a 

o 

a 

z 

O  LU 

X 

a 

a.  o 

ID 

in 

X 

a  • 

z  z 

0 

2 

CN  * 

a 

Z  a 

LU 

tn 

in 

CJ 

CN 

a 

x  a 

a 

*— i 

•  0 

< 

Q 

LU 

z 

0 

Z  a. 

o 

H  O 

1  D 

a 

0  a 

a 

II 

a 

Z 

X 

X 

< 

z 

uj  2 

H 

o 

a 

D 

a 

in 

in 

Q 

a  hi 

♦— 

0  3 

hh 

K 

u 

in 

< 

X  < 

h- 

UJ 

a  a 

<  -i 

a 

o 

a 

£1 

Z  2 

Z  -J 

CD 

in 

1-  ex 

o 

a 

o 

O  a 

3 

o 

in  - 

o 

o 

CJ 

*-«  3 

0 

cn 

*  < 

< 

CO 

0 

O 

b- 

o 

a  = 

2  K  H 

a 

a 

*  X 

a 

CJ 

00  CJ 

Z 

LU 

CJ  3 

< 

O  ex 

0 

3 

a 

O 

3  Q 

2 

< 

0  a 

o  a 

o 

z 

C/1  CJ 

*— i 

K 

CJ  b- 

ll 

•  < 

0 

a 

H  O 

+  0  O 

♦— t 

> 

a  2 

TT  HI 

0 

o  < 

in 

< 

Z  CJ 

O  LU 

O 

Z 

2 

3  in 

Z  2 

O 

a  o 

a 

•f 

ex 

in 

-J 

-  < 

t— t 

uj  2 

H- 

• 

HH 

0 

1-  • 

a 

a 

cn  o 

cn  - 

a 

h- 

CJ  < 

o 

3 

O 

O 

* 

■*“ 

•  HH 

6 

o 

z  o 

3  O 

a 

2  O 

\ 

<[ 

o 

ex 

2 

O 

O 

i 

f-  -J 

UJ 

o  o 

o 

0  O 

in 

0 

2  > 

a  O 

* 

CJ 

LO  2 

0 

3 

cn  * 

CN 

* 

*— * 

UJ 

* 

3 

UJ 

«- 

a  = 

2  r-  * 

a 

a 

*—» 

2 

•-»  a 

H-t  LO 

w  UJ 

r-> 

z 

CJ 

01  w 

0) 

> — * 

* 

3 

w ' 

0  UJ 

K 

01 

in 

01  W 

3 

3 

X 

< 

a 

a  cn 

3 

l/) 

CJ 

a 

h- 

CN 

2 

b- 

Z  X 

3 

X 

H  2 

3  a 

l»  K 

Z 

Z 

a 

C3  2 

3 

a  2 

II 

u  * 

•"* 

< 

a  < 

b- 

< 

*—♦ 

< 

a 

0 

3 

b- 

< 

3  X 

K  < 

>— ( 

a 

a  a 

a 

«3  HH 

CJ  < 

Z  2  Z 

2 

CJ 

h- 

2 

U. 

z 

*— i 

Z  2  ll 

Q  a 

K  Z  2 

1— 

a 

a 

in 

ex  in 

a  2  ^  a 

< 

Z  lu 

X 

CO 

*-<  ex 

HH 

ex 

LU  2 

ex 

~  o 

w 

N— ' 

HH 

a 

< 

3  h  ex 

Z 

Z 

< 

< 

o  a 

a  hh 

a  2 

CJ 

ex 

CJ 

CJ 

ex  o 

ex 

o 

a 

O 

o 

z 

Ll. 

LU 

ex 

o  < 

Z 

0  a  o 

o 

o 

X 

a  a 

<  a 

O  D 

in 

z 

< 

z 

z 

a  u_ 

CL 

LU 

in 

CJ 

Lu 

< 

a 

a  2 

*— i 

Z  a  a 

0 

0 

a 

in  a 

0  a 

a  0 

0 

O 

o  o  o 

c 

O  O 

o  o  o 

O  O 

<. 

O 

o  o  o 

o 

O  in 

o  o  o 

O  O 

a 

cn 

CN 

O' 

in 

in 

^  CN 

h«- 

cn 

0 

in  cn 

H* 

tn 

0> 

cn 

01 

01 

Ol 

C 

C 

C 

CJ 

c 

c 

CJ 

c 

c 

c 

0  0  0  0 

C 

C 

0 

C 

C 

0 

■r* 

I* 


< 

a 

0 

o 


o 

CN 


o 

cn 


in 

cn 


O 

TT 


in 


o 

in 


in 

in 


B-4 


SPLINE  FITS  SIX  SPECTRAL  POINTS  WITH  TWO  CUBIC 
SPLINE  CURVES:  ONE  FOR  THE  REAL  PART  AND  ONE  FOR 
THE  IMAGINARY  PART.  THE  CURVATURES  ARE  STORED  IN 
THE  ARRAY,  G. 


REPORT  NO.  NADC-82253-30 


LU 

O 

< 

CL 


o 


tn 

in 


r- 

O 


<T> 

o 


CL 

in 


CM 

ID 

O 


o 

in 

O 


a 

in 


cc 

< 

CL 


CN 

(I 

o 


# 

DC 

a 

"3 

♦ 

X T 

< 

Z3 

l 

* 

2 

li 

*• 

CM 

i— i 

£ 

N 

+• 

a 

O 

O 

* 

< 

< 

H 

1 

* 

CL 

£ 

K 

“3 

* 

•  ^ 

►— < 

o 

w 

#— s 

* 

~  CD 

O 

OQ 

> 

CM 

* 

O  w 

DC 

-w' 

* 

in  > 

O 

> 

£ 

LU 

+ 

O  in 

* 

o  • 

+ 

O 

•  < 

* 

CN 

+ 

DC 

II 

o\ 

* 

w  ID 

►— < 

Lu 

O 

”3 

* 

<  w 

II 

* 

f— \ 

CD 

■w 

II  CM 

* 

k  in 

o 

CM 

in 

tn 

"3 

> 

* 

<  LU 

1— « 

o 

w 

2 

w 

* 

tn 

* 

O  • 

w 

+ 

> 

O 

in  o 

o 

a  lu 

# 

"<T 

* 

*— i 

LU 

H 

< 

Vh 

« 

^  in 

CD 

-J 

<r- 

o 

H- 

\  CM 

'  -  II 

* 

Tt  W 

CN 

< 

•— i 

II 

< 

O 

o 

ID 

* 

w  m 

i 

LU 

CD 

w 

CM 

ZD 

o 

1 

* 

—  < 

DC 

* 

< 

< 

O 

CD  O 

* 

O 

* 

II 

K 

i 

LU 

/— s 

^-s 

t— i 

* 

<  2 

* 

o 

h- 

ll 

< 

< 

T- 

1 

o 

n 

* 

CL  O 

“3 

►—< 

o 

•— « 

Q 

'—V 

LU 

+ 

CM 

CM 

* 

>— < 

★ 

< 

ID 

O 

ii 

6 

O  “3 

•*— 

/— s.  w 

* 

2  in 

CN  O 

DC 

O 

II 

ZD 

LU 

w 

O  CD 

* 

O  2 

in 

H 

■*— 

in 

> 

Q 

TT 

> 

o 

>— i 

-* 

£  lu 

11 

*“ 

X 

LU 

"3 

"3 

H 

>  ll 

* 

£  £ 

LU 

ii 

II 

DC 

W 

II 

II 

* 

O  f-i 

o 

o 

w 

tn 

Lu 

O 

w  O 

* 

U  o 

*-h 

D 

Q 

> 

"3  lu 

< 

*— i 

< 

Lu 

CD 

a  cd 

o 

cd 


O 
CD  ' 


O 


O  LU  2 
CD  OC  UJ 


o 

in 


*  *  *  * 


t- 

z: 

o 

DC 

QQ 

D 

in 


in 

CM 


o 

CD 


in 

CD 


B-5 


COMPUTE  CURVATURES  FROM  TOP  DOWN 


REPORT  NO.  NADC-82253-30 


UJ 

CJ 

< 

0. 


■O 


in 

in 


O 

N 

01 

o 

N 

CN 

oo 


CD 

CN 

TT 

+ 

10 

Z 


Li¬ 


lt 


;  cl 
fc'0 


r- 

N 

x 


LU 

CO 

< 

I 

Q. 

UJ 

Z 


UJ 

Z 


2 

< 

a 

UJ 

00 

< 

X 

a 


X 

o 

dc 

m 

3 

to 


00  LU 


z 

UJ 

2 

_ 

a 

Z 

on 

o 

X 

i— i 

»— « 

in 

X 

X 

l-H 

a 

» 

00 

LU 

in 

UJ 

X 

O 

X 

Q 

X 

o 

00 

x 

cj 

o 

X 

o 

CL 

* 

tH 

X 

LU 

Q 

U 

DC 

LU 

Q 

LU 

X 

X 

DC 

X 

< 

CJ 

X 

O 

X 

cj 

X 

X 

X 

ac 

X 

dc 

X 

z 

Z 

CQ 

X 

X 

CJ 

a 

< 

2 

l-H 

CL 

LL 

oo 

X 

Z 

X 

on 

O 

x 

l— ( 

a 

l-H 

Z 

t-H 

< 

< 

u_ 

o 

LU 

X 

QG 

o 

» 

a 

2 

X 

o 

X 

X 

O 

X 

X 

X 

X 

UJ 

00 

DL 

z 

X 

x 

Z 

00 

LU 

X  LU 

O 

< 

a 

►H 

00 

* 

z 

X 

X 

X 

z 

DC 

X 

o 

* 

* 

*— ( 

o 

o 

l-H 

O 

X 

DC 

* 

* 

X 

< 

LU 

z 

a 

X 

o 

* 

* 

X 

> 

X 

LU 

X 

* 

* 

o 

z 

00 

O  X 

< 

oo 

oo 

* 

-» 

DC 

O 

»— t 

LU 

x 

X 

X 

t-H 

z 

* 

* 

t— t 

X 

X 

z 

o 

« 

* 

UJ 

00 

00 

X 

00 

- 

< 

X 

on 

t-H 

* 

* 

> 

on 

X 

MH 

HH 

o 

t-H 

X 

* 

* 

M 

UJ 

Z 

O 

X 

(J 

cj 

X 

CJ 

+ 

* 

X 

a 

►h 

z 

LU 

X 

< 

X 

* 

* 

X 

o 

o 

HM 

_l 

t-H 

a 

X 

DC 

* 

in 

* 

CJ 

o 

CL 

00 

cj 

on 

X 

X 

t-H 

* 

CN 

* 

LU 

cr 

X 

z 

cj 

< 

o 

* 

o 

* 

X 

a 

X 

< 

o 

o 

z 

X 

* 

* 

UJ 

z  >- 

CL 

Z 

DC 

/— S 

* 

* 

UJ 

LU 

CQ 

oo 

X 

X 

* 

C J 

* 

UJ 

00 

cj 

i-h 

IN 

• 

X 

00 

* 

X 

< 

* 

X 

< 

< 

X 

X 

> 

t-H 

* 

X 

X 

* 

X 

X  X 

X 

Z  in 

00 

l-H 

3 

* 

u 

X 

* 

CL 

D  LU 

< 

X 

< 

X 

* 

* 

z 

* 

00 

< 

LU  X 

X 

X 

< 

CJ 

* 

\ 

* 

►-H 

_l 

CJ 

O 

x 

DC 

a 

cj 

o 

* 

\ 

X 

* 

< 

O  Z 

CJ 

X 

X 

* 

X 

00 

* 

LU 

a 

z 

< 

LU 

DC 

X 

X 

Z  o 

* 

X 

X 

* 

in 

x 

l-H 

z 

X 

a 

X 

1 

* 

(J 

X 

* 

< 

u 

X 

Q 

l-H 

X 

X 

X 

DC 

* 

\ 

\ 

* 

X 

UJ 

on  uj 

00 

< 

o 

X  x 

* 

* 

a 

CL 

UJ 

Q 

O 

X 

01 

3 

X 

X 

* 

z 

z 

* 

on 

X 

X 

CJ 

DC 

a 

z 

* 

o 

o 

* 

_J 

cj  z 

X 

X 

X 

* 

2 

2 

* 

LU 

DC 

cj 

UJ 

< 

X 

< 

o 

* 

2 

2 

* 

X 

o  z 

X 

00  X 

o 

X 

o 

* 

o 

o 

* 

i- 

x 

HH 

x 

t-H 

X 

X 

X 

* 

u 

u 

*  * 


X 

on 

D. 

Q 

X 

00 

z 

<***" 

Z 

oo 

< 

CN 

t-H 

t-H  \ 

X 

CN 

z 

* 

o 

t-H 

CJ 

> 

z 

o 

10 

CL 

00  x 

X 

- 

t-H 

►h  a 

DC 

CN 

< 3 

X 

cj 

X 

X  \ 

X 

X 

t-H 

< 

- 

O 

X 

Tf 

oo 

DC 

CJ 

/“s 

o 

z  o 

O 

o 

o 

in 

x  2 

c 

o 

z 

X 

on 

in 

DC 

»-*  < 

CJ 

< 

X 

CJ 

O 

D 

X 

X 

>- 

X 

CN 

X 

X 

• 

. 

-- 

X 

DC 

DD 

* 

CJ 

•  > 

CN 

00 

> 

o 

U 

X 

t-H 

X 

- 

o 

X 

z 

X 

H— v 

> 

t-H  X 

o 

X 

X 

CJ 

3 

a 

CN 

X  ~ 

DC 

X 

X 

o 

00 

+ 

z 

CL  00 

X 

O 

X 

X 

a 

• 

in 

X 

O 

—  o 

x  Z 

CJ 

< 

X 

* 

t-H 

00  x 

- 

o 

o 

< 

DC 

cj 

Z 

CN 

X 

o 

X 

□ 

X 

z 

t— < 

t-H 

CJ  X 

in 

O  Q 

DC 

HH 

CN 

o 

< 

in  \  co 

- 

oo 

X 

a 

X 

Q 

Q. 

X 

o 

o 

DC 

CL 

z 

u 

< 

CN 

< 

0-  \  X 

in 

X 

X 

< 

< 

a 

w 

X 

O 

\  00 

X 

l-H 

X 

X 

2 

DC 

* 

< 

< 

X 

<  X 

X 

* 

z 

X 

oc 

H-K 

X 

X 

DC 

ii 

I\Z 

DC 

H—S 

X 

o 

X 

2 

o 

< 

X 

X 

CJ 

CN 

HH 

CJ 

X 

DC 

in 

O 

CJ 

CN 

11  ^ 

X 

> 

CN 

o 

X 

X 

X 

< 

O 

CN  o 

X 

>- 

X 

t-H 

o 

- 

> 

o 

Cl 

in 

z 

o 

O 

X 

>  CN 

It 

a 

o 

m  X 

CN 

X 

-  'W 

in 

•  CL 

o 

X 

CN 

Z 

o 

O 

CN 

Q 

CN  < 

ifi 

X 

<  W 

01 

. 

X 

o 

X 

a 

X 

2 

X 

X  X 

> 

X 

(J 

in 

CN 

T— 

X 

X 

o 

> 

in 

-  < 

o 

X 

o  z 

> 

>- 

a 

a 

X 

X 

X 

—  o 

a 

••  CJ 

< 

* 

- 

> 

o 

z 

00 

DC 

CN 

> 

DC 

tt  h- 

a 

CN 

X 

•*- 

▼— 

o 

< 

X 

X 

X 

2 

> 

-  ^  H- 

t-H 

CN  •+ 

DC 

>  Z  X 

> 

X 

H— 

o 

X 

K 

X 

X 

rr  w 

< 

O  X 

CL 

-  x 

- 

X  o 

a 

in 

X 

X 

X 

/— v 

z 

X  w  < 

D. 

* 

CN  X 

in 

+ 

X 

o 

X 

CJ 

X 

X 

o 

-  o 

H— 

‘  CN 

X 

X  X 

in 

. 

X 

li 

o 

X 

X 

X 

\  Q 

X 

v 

o 

X 

CN  O 

o 

X  X 

X 

X 

z 

Q. 

Cl 

in 

X  <  z 

< 

o 

It  < 

o 

o  - 

_ j 

X 

o 

X 

X 

X 

o 

t-H 

\CLOOX 

< 

X  X 

z  x 

* 

* 

H- 

X 

u 

X 

X 

< 

X 

t-H 

< 

X 

< 

in 

<  > 

X 

T- 

oo 

w 

+ 

u 

Z 

X 

X 

CL 

a 

X 

CJ 

Z  Z  00 

It 

Q 

o  o 

t-H 

t-H 

X 

X 

X 

X 

t-H 

o 

X 

X 

00 

CJ 

o 

< 

o  o  z 

X 

o 

X 

X 

< 

in 

X 

X 

X 

X 

2  2  x 

X 

H 

in 

ID  li 

in 

X  t-t 

X 

IS 

z 

2 

CJ 

(J 

X 

o 

< 

CL 

X 

X 

X 

X 

2  2  2 

X 

X 

CJ 

•  in 

oo 

l-H 

DC 

< 

t-H 

DC 

CJ 

2 

X 

X 

X 

O  O  ^ 

CJ 

X 

O  CN 

^  O 

X 

in 

DC 

O 

X 

X 

X 

a 

O 

< 

< 

< 

UOQZ 

X 

O  X 

X  CL 

2 

u 

CL 

X 

X 

l-H 

X 

X 

CJ 

o 

CJ 

u 

<  <  X  O 
(J  O  X  CJ 


* 

* 


*  +  *  * 


******  * 


ouu 


o 

o 

01 

uoo  uuuouu  o 


c 

o 

ro 

ouuouu  o  cj  cj 


X 

o 

dc 

DQ 

X 

00 


in 


O 


O 

CN 


in 

CN 


O 

ro 


in 

ro 


O 


in 

^r 


C  in 

in  in 


B-6 


CHECK  FOR  CROSSING  OF  THE  NEGATIVE  REAL  AXIS. 


SUBROUTINE  PHASE  74/74  OPT=1  FTN  4.6+428  82/09/07.  15.51.04  PAGE 


REPORT  NO.  NADC-82253-30 


CO 

h- 

Z 

X 

NH 

< 

o 

Cl 

< 

X 

LU 

H 

a 

s 

Ui 

+ 

> 

X 

o 

< 

K 

c 

LU 

CL 

z 

X 

LU 

CO 

I 

C3 

H 

Z 

U. 

CO 

O 

in 

>— « 

o 

C3 

CN 

oi 

z 

> 

u 

CN 

t— < 

• 

f- 

> 

m 

CN 

LL 

X 

in 

X 

O 

o 

II 

o 

• 

z 

01 

■»— 

01 

•*- 

u 

> 

UJ 

II 

> 

• 

CD 

a 

■»- 

£ 

o 

X 

X 

CN 

LL 

'w' 

Z 

+ 

h-  UJ 

X 

LU 

* 

z>  x 

UJ 

-ft 

CN  X 

o 

u  z 

01 

CN 

X  z  z 

LU 

1— 1 

O 

w 

*-t  01 

I 

_J  h- 

H 

o 

II  h-  X 

o 

-J  z 

in 

UJ 

z  k  a 

<  o 

CL 

O  UJ  Z 

CJ  u 

CO 

x  a  oi  lu 

o  o 

o 

o  o 

o 

rr  in 

ID 

U  UU  UUU 


o  in  o 

id  id  t" 


B-7 


CUT  TESTS  A  PAIR  OF  SPECTRAL  POINTS  TO  DETERMINE 
IF  THE  NEGATIVE  REAL  AXIS  HAS  BEEN  CROSSED  AND,  IF  SO, 
WHAT  THE  DIRECTION  OF  CROSSING  IS.  THE  CASE  OF 
Y 1  =  Y 2  -  O  HAS  BEEN  IGNORED. 


REPORT  NO.  NADC-82253-30 


U 

< 

CL 


o 

in 

in 


* 


o 

Oi 

o 

CN 

00 


b 


00 

CN 

TT 


X 

U 


x 

u 

z 


x 

u 


z 

o 

E 

E 

Q 

CJ 


a 

z 

h- 

O 

X 

o 

»— i 

l/) 

t— 1 

* 

LU 

a 

LU 

i- 

DC 

Z 

X 

O 

*— * 

CJ 

LU 

z 

a 

to 

C 

01 

H* 

a 

Z 

UJ 

1- 

t-H 

►- 

-J 

o 

z 

z 

CJ 

a 

cc 

a 

LU 

X 

X 

►“ 

LU 

O' 

t- 

h* 

LU 

CQ 

LU 

LU 

LU 

Q 

X 

QC 

a. 

CN 

1- 

t— 

O 

l>0 

1— 1 

/— * 

1 

z 

x 

LU 

o 

O 

4-4 

X 

L— 

to 

U- 

o 

o 

X 

to 

LU 

HH 

o  ■ 

CJ 

o 

— 1 

►- 

O  t- 

z 

tx 

o 

LU 

o 

-  a 

CJ 

z 

-J 

z 

II 

< 

CQ 

CN 

a  o  cn 

X 

*— i 

X 

X  H  > 

H 

CJ 

Q 

to 

X 

z 

LU 

to 

O' 

iu  o  a 

CJ 

h- 

< 

a 

O 

O 

a  o  o 

Z  X  a: 

x 

CL 

• 

+ 

U 

CQ 

_j 

o 

o 

z 

cj 

to 

o  o  • 

CN 

• 

z 

<— t 

LU 

o 

•  •  O  X  >  O 

1—4 

z 

CN 

o  o  • 

CJ 

o 

o 

< 

h- 

> 

■  •  ►—  z 

l~  O 

* 

LU 

z 

— 1 

* 

hh  u 

— 1 

C T> 

w 

.CJ 

1— 1 

a 

O  -J  • 

II 

K 

Z 

to 

> 

.  .  1— 

*— 

H* 

< 

to 

u. 

X 

CL  D.  > 

1“ 

> 

z  s 

to 

o 

-j 

w 

(1 

W  W  > _ ' 

X 

w 

1-4 

a 

cn 

< 

LL 

LLlLlL 

CJ 

U_ 

a 

o 

u  X 

>— t 

a 

►— <  4— 4  (— 1 

z 

*— ( 

Cl 

U_ 

K  D 
lu  Z 


O 

O 


uuuuu 


o 

o 

o 

a> 


3 

o 

cz 

CD 

x 

to 


in 


O 


O 

CN 


tn 

CN 


f 


B-8 


SUBROUTINE  CPA(P,J) 


REPORT  NO.  NADC-82253-30 


Llj 

(D 

< 

Q_ 


* 


Tt 

O 


in 

in 


O 

<n 

O 

CM 

00 


00 

CN 

rr 

+ 

10 


< 

CL 

u 

LU 

2 


< 

2 


LU 

2 


o 

< 

2 

DC 

DC 

O 

CL 

H 

< 

PH 

Cl 

2 

CL 

h- 

< 

LU 

CJ 

2 

2 

UJ 

UJ 

10 

O 

tn 

I 

LU 

«— i 

PH 

K 

in 

h* 

CQ 

<* 

U. 

LU 

CJ 

> 

o 

2 

o 

00 

p- < 

_J 

UJ 

-J 

~D 

T- 

2 

3 

CL 

• 

tn 

+ 

* 

-J 

in 

* 

o 

PH 

3 

o 

* 

< 

* 

* 

> 

LU 

* 

^p* 

a 

CL 

* 

2 

* 

w 

H- 

2 

O 

a 

* 

DC 

K 

* 

K 

in 

< 

* 

LU 

* 

in 

*— i 

O 

UJ 

* 

b- 

Lu 

* 

•— t 

a 

3 

t~ 

Q 

* 

LU 

O 

# 

a 

* 

2 

* 

LU 

in 

UJ 

+ 

* 

< 

2 

* 

ii 

o 

*- 

U 

* 

DC 

cj 

* 

2 

2 

CL 

* 

< 

< 

* 

CN 

K 

PH 

in 

< 

* 

CL 

o 

* 

Q 

2 

o 

h- 

o 

II 

* 

cr 

* 

UJ 

CL 

o 

in 

o 

* 

LU 

a 

* 

2 

PH 

CL 

* 

2 

CL 

# 

* 

3 

2 

n 

o 

* 

h- 

< 

* 

O 

LU 

a  o 

* 

* 

CL 

LU 

CL 

h 

-  h- 

* 

i n 

LL. 

* 

:  < 

3 

in 

CN 

* 

LU 

o 

* 

3 

H- 

UJ 

Q  O 

o 

* 

H 

* 

• 

H 

UJ 

v*  O 

•  CJ 

« 

< 

h- 

* 

O 

in 

CQ 

_i 

3 

3 

CN 

* 

CJ 

2 

# 

ce 

- 

< 

o  ~ 

a  \ 

* 

O 

p- • 

* 

o 

i— « 

LU 

in  *- 

2 

•  CM 

CL 

CL 

CL 

* 

O 

* 

w 

LU 

> 

CM  II 

in 

o  a 

LU 

_j 

* 

a 

* 

h- 

DC 

•  a 

o  • 

K  O  H 

o 

UJ 

* 

< 

* 

in 

a 

3 

O  m 

UJ 

O  H 

in  in 

in 

O 

« 

CL 

i— 

+ 

PH 

2 

CJ 

K 

CD  _l 

PH  H~ 

PH 

1 

* 

cj 

tn 

* 

a 

< 

ii  o 

< 

Q 

a 

n 

* 

UJ 

* 

w* 

CO 

o 

CJ 

h* 

O 

a 

* 

in 

* 

II 

2 

a  cm 

O 

2  Q 

It  K 

n 

CL 

* 

o 

# 

a 

o 

-j 

-1 

PH  W 

n 

_l 

* 

-J 

* 

«*- 

-J 

LU  O 

DC  LU 

■*-  O 

CN 

UJ 

* 

o 

* 

o 

< 

a  d 

CL  ph 

a  o  o 

a 

Q 

* 

* 

« 

-# 

o 

O 

* 

+ 

o 

in 

* 

* 

* 

+ 

+ 

*  * 

*  * 

* 

C 

C 

c 

c 

C 

c 

c 

CJ 

2 
O 
DC 
CQ 
ZD 
i n 


in 


O  in 

CM  CM 


B-9 


200  CONTINUE 

9000  FORMAT ( *  D1  D2  P  *3F20.4) 
RETURN 
END 


DIST  FINDS  THE  (DI STANCE ) *+2  TO  THE  ORIGIN  OF 
POINT  ON  A  SINGLE  SPLINE  SEGMENT. 


REPORT  NO.  NADC-82253-30 


a 

< 

CL 


O 


in 

in 


O 


CN 

CO 


CO 

CN 

-ft 

ID 


LO 

O 


Z 

o 


CN 

ID 

o 

O 

is 

*■" 

CD 

CD 

o 

ft- 

h- 

in 

ft 

* 

O 

CO 

^  CD 

/-V 

CN 

f- 

■—  h- 

CN 

W 

w 

-ft  w 

-ft 

o 

* 

"0  ft 

"D 

UJ 

* 

ft 

a 

CN  CN 

CN 

in 

• 

w 

- 

<  -O’ 

< 

W 

ft-  w 

►- 

Tt 

o 

<  a 

< 

w 

ft 

o  * 

O 

CN 

00 

*  CO 

« 

a 

K 

CO  ft- 

CD 

< 

h- 

ft- 

a 

-ft 

+ 

* 

+ 

-ft 

O  CD 

o 

^  o 

in  h* 

o 

1  ■»— 

CN  l 

1 

“D  » 

“D 

■w' 

< 

ft  < 

ft 

<  O  H 

CN  ft- 

CN 

H  • 

-ft 

N— '  ft 

•w^ 

<  — 

< 

<  < 

< 

Q 

H 

ft-  ft- 

1- 

> 

*  ii 

-w 

<  w 

< 

ft 

ft 

Q  ft 

a 

> 

Tf  < 

S—' 

ft 

-* 

w  H 

■»— 

<  CN 

< 

-ft 

-- 

- 

ft-  - 

H 

a 

CO 

CO 

X 

<  ^ 

w 

-ft  w 

+ 

-*■ 

a 

CJ 

o 

X 

ft 

ft 

z  a 

< 

< 

1! 

o 

1— 

ft- 

S  II 

ft- 

s 

11 

II 

in 

O  CD 

►— * 

a  ►— 

x 

> 

Q 

** 

*  *  >— 


a 
z 

o  a  lu 


u 

z 

3 


B-10 


SUBROUTINE  DFT  74/74  OPT=1  FTN  4.6+428  82/09/07.  15.51.04  PAGE 


REPORT  NO.  NADC-82253-30 


2 

cl 

o 

u. 

CJ 

Ul 

z 

Z 

< 

< 

H 

a 

h- 

> 

m 

cl 

LU  ■ 

in 

Q 

*-*  a 

O 

UJ 

LU 

or 

■*- 

«_ i 

★ 

X  D 

* 

HH 

o : 

< 

-* 

O  Z 

¥■ 

• 

UJ 

f— 

-* 

LU  < 

+ 

CJ 

in 

o 

♦ 

* 

z 

a: 

* 

UJ  Z 

* 

o 

< 

UJ 

2 

* 

H- 

* 

CM 

H 

2 

in 

UJ 

* 

# 

UJ  > 

* 

U_ 

w 

LU 

_j 

< 

* 

CL  CO 

* 

* 

Z 

K 

C5 

a 

* 

u 

* 

CL 

Z 

• 

♦ 

CO  Q 

* 

in 

D 

UJ 

< 

* 

*-•  UJ 

* 

UJ 

CO 

K 

* 

D  hh 

* 

H 

o 

H  O 

f- 

Ll 

* 

* 

u. 

* 

ll 

Ll 

t- 

o 

or 

* 

UJ  1-. 

* 

Q 

c  o 

in  o 

t_ i 

u. 

Z  UJ 

H 

-* 

X  U 

* 

w 

Q 

X 

O 

in  z 

in 

* 

* 

H*  UJ 

* 

00 

00  II 

n 

to 

o  in 

o 

Z 

* 

a 

* 

O 

TT 

Z 

z 

*  * 

(J 

• 

* 

in  oo 

* 

ro 

O  uj 

«  2 

z 

»— * 

0. 

* 

UJ 

* 

Z 

CN  Z 

z  z 

o 

s— ' 

* 

►—  Z 

* 

•jS 

N  in 

Q 

z 

2 

z 

1- 

* 

x  o 

* 

Z 

< 

a 

ID 

CL 

<  < 

< 

u. 

* 

a.  •“« 

* 

n 

► 

CL 

* 

0 

in 

UJ 

o  a 

Q 

* 

2  H- 

* 

< 

00 

CL 

CL 

*-h  m 

Z  W 

i- 

CN 

H- 

z 

* 

o  < 

* 

o 

CN 

* 

+ 

a 

< 

LL 

. 

-1-  1 

UJ 

* 

<J  CJ 

* 

o 

o 

f- 

Q 

y— 

K 

CO 

Z 

* 

o 

* 

z 

CD 

o 

V 

s  o 

w  O 

II 

X 

UJ  2 

i— i 

* 

f-  _J 

* 

o 

* 

■*“ 

h-  • 

in  . 

2 

z 

UJ 

cr  •-< 

UJ 

h* 

+ 

u. 

* 

1— 1 

II 

O) 

1 

T— 

o  o 

CL 

Z 

o  a 

X 

* 

D  uj 

* 

CO 

H 

z 

II 

u 

□ 

o 

< 

O 

* 

X 

* 

z 

t— ( 

t— 

< 

w 

II 

If 

ll 

o 

Q 

ii  it 

o 

a 

* 

t- 

* 

UJ 

CL  Z 

2 

cj 

II 

CL 

in 

Q 

Cl 

m 

+ 

* 

2 

O 

a 

II 

z  z 

UJ 

LU 

< 

uj  2 

D 

X 

* 

h- 

* 

HH 

o: 

o 

<  in 

cj  a. 

CL 

o 

C£  hh 

CO 

* 

< 

* 

Q 

H 

a 

ll 

< 

1-  o 

a  a 

a 

Q  Q 

* 

* 

* 

* 

O 

*■ 

* 

o 

* 

•* 

1— 

* 

* 

CD 

* 

*  *  *  * 

o 

V 

c 

c 

CJ 

c 

c 

CJ 

O 

CO  in 
Q  O 
*  * 

UJ  z 
Z  co 

CO  U 
I  + 

cj  cj 

D  O 
*  * 

Z  uj  hi 
co  Z  X  lu 
CJ  CO  (J  Z 

z 

It  <1  II  >-l 

h- 

uj  uj  Z  2 
x  Z  in  o 
<J  in  o  CJ 

O 
O 
in 

O  o 


CLCLHI—  aUQCUJ 


CJ 


o 

o 

o 

cn 


o 

cn 


in 

CN 


O  in 

co  ro 


o 


B-ll 


COMPENSATE  FOR  SHIFT  IN  TIME  SERIES 


ACUT  TESTS  A  PAIR  OF  SPECTRAL  POINTS  TO  DETERMINE 
IF  THE  NEGATIVE  REAL  AXIS  HAS  BEEN  CROSSED  AND,  IF  SO, 


REPORT  NO.  NADC-82253-30 


LU 

O 

< 

a 


TT 

o 


in 

in 


r- 

C 

N 

O) 

o 

\ 

CN 

00 


s 

o 


oi 

LU 


3 

CJ 

< 


‘  LU  LU 

u- 

o 

O  ON 

o 

a 

z 

ON 

UJ 

* 

(/)  UJ 

in 

<  Z 

ON 

CJ 

CJ  ON 

CJ 

K 

z 

UJ  X 

X 

X  o 

u 

k  o: 

z 

o 

in 

a 

UJ 

•  ON 

LU 

* 

H* 

in  x 

t- 

* 

z 

z 

CJ 

ON  K 

z 

* 

01 

01 

LU 

X 

* 

X 

X 

h- 

CJ 

o 

* 

h- 

H* 

U> 

z  • 

u 

+ 

LU 

UJ 

a 

ON  O 

* 

01 

01 

in  uj 

* 

CJ 

<n  o£ 

< 

* 

/— s 

CN  Z 

o  o 

CO 

* 

o 

o 

ON 

01  2 

o 

* 

i  in 

CJ  O 

_J 

* 

o 

o 

in 

ON 

o 

* 

o  • 

CJ  o 

LU 

* 

h- 

O  t— 

CJ  01 

o  z 

UJ 

* 

CJ 

—  CJ 

z  o 

UJ 

X 

* 

z 

Z  UJ 

K 

* 

CN 

01  O  CN 

II  X 

o  m 

( 

* 

X 

X  K  > 

CJ 

ON 

►“  1 

;  * 

h* 

CJ  Z 

t-  in 

<  1 

1  * 

01 

LU  O  01 

CJ  < 

cj  < 

X 

* 

o 

01  CJ  o 

Z  CJ  01 

UJ  X 

♦- 

* 

CJ  CD 

m 

* 

o 

o 

~  Z 

ON  O 

z 

* 

o  o  •  + 

CN  ■> 

o 

ON 

* 

o 

•  ■  o 

>  o 

UJ  >  O 
I  CN  _J  LU 
t-  >  Z  l“H 
O  U_ 
h  II  t-M 


< 

X  *•“  X 
3  >  CJ 


*  ■  CN  o  O  •  CJ  •  O 

*  1“  >-  ■  •  H  CJ  J-  O  * 

*  CJ*  hh(J2  JOIw 

*  •  —  O  -I  •  •  t-  2 

*  ^  >  *  •  —  it  ■*—♦—<  QC 

*><  Q.  Q.  >  >  Z  2  X 

*  w  n  wwwUwt-iat-Q 

*  u.  L.L.iLuu.aoiu2 

*  mQ.mi-(h2wQ.ILQ!UJ 


O 

o 


*  *  *  *  * 


o 

o 

o 

cn 


X 

o 

a 

co 

x 

10 


in 


O 


O 

CN 


3-12 


REPORT  NO.  NADC-82253-30 


No.  of 
Copi es 

Defence  Research  Centre  1 

Salisbury,  Adelaide,  SA,  Australia 

Defence  Scientific  Establishment  1 

HMNZ  Dockyard 

Auckland  9,  New  Zealand  (Dr.  K.M.  Guthrie) 

DTIC  12 


> 


REPORT  NO.  NADC-82253-30 


No.  of 
Copies 


Applied  Physics  Laboratory  1 

The  John  Hopkins  University 
John  Hopkins  Road 
Laurel,  MD  20707 

Lamont-Doherty  Geological  Observatory  1. 

Palisades,  NY  10964 

Woods  Hole  Oceanographic  Institution  1 

Woods  Hole,  MA  02543 

Scripps  Institute  of  Oceanography  1 

Marine  Physics  Laboratory 
San  Diego,  CA  92152 

Auditory  Research  Laboratory  1 

Northwestern  University 
2299  Sheridan  Road 

Evanston,  IL  60201  (Dr.  Frederic  Wight*  r 

Mayo  Medical  School  1 

Rochester,  MN  55901  (Dr.  James  Greenleaf) 

Canadian  Superior  Oil  Ltd.  '  1 


Three  Calgary  Place 
355-4th  Avenue  SW 

Calgary,  Alberta  T2P  0J1  (Dr.  R.  Hinds) 


Defence  Research  Establishment,  Pacific  1 

Forces  Mail  Office 
Victoria,  BC  Canada  VOS  1B0 

Defence  Research  Establishment,  Atlantic  1 

Forces  Mail  Office 
Halifax,  NS  Canada 

Admiralty  Underwater  Weapons  Establishment  1 

Portland,  Dorset,  England,  UK 

Admiralty  Research  Laboratory  1 

Teddington,  England,  UK 

Royal  Aircraft  Establishment  1 

Farnborough,  Hants,  England,  UK 

Royal  Australian  Naval  Research  Laboratory  1 

Edge hi  11  Road 

Sydney,  NSW,  Australia 


REPORT  NO.  NADC-82253-30 


DISTRIBUTION  LIST 
TASK  AREA  NO.  ZR  01108 


No.  of 
Copies 


CNM 


(1  for-ASW-01 ) 

(1  for-ASW-10) 

(1  for-ASW-13) 

(3  for-MAT-0721 ) 


9 


(3  for-MAT-0812) 
NAVAIRSYSCOM  (00D4) 


6 


(2  for  retention) 

(1  for-AIR-340) 

(1  for-AIR-950D) 

(1  for-AIR-54902B)' 

(1  for=PMA-264) 

0NR-102-0S  1 

NAVELEXSYSCOM-PME-124  1 

Naval  Postgraduate  School  (Prof.  0.  Heinz)  1 

NOSC-714  1 

N0RDA-320  J, 

NRL-8100  1 

NUSC,  New  London  1 

SACLANT  ASW  Research  Centre  1 

Viale  San  Bartolomeo  400,  1-19026 
La  Spezia,  Italy 

Applied  Research  Laboratory  3 

The  Pennsylvania  State  University 
P.0.  Box  30 

State  College,  PA  16801  (Dr.  S.T.  McDaniel;  Dr.  L.  Sibul;  library) 


